// ACA.cpp : Defines the entry point for the console application.
//
extern "C"
{
#include "induct.h"
}
#include <cmath>


/*Here the intersection between a line and a plane is calculated*/
void xyz_point(double x_y_z[3], double grad[3], double xyz_obs[3], double x, double y, double z)
{
	double al = pow(grad[0], 2);
	double bm = pow(grad[1], 2);
	double cn = pow(grad[2], 2);
	if (grad[0] != 0)
	{
		double bm_plus_cn = bm + cn;
		x_y_z[0] = (x * al + xyz_obs[0] * bm_plus_cn - grad[1] * grad[0] * (xyz_obs[1] - y) - grad[2] * grad[0] * (xyz_obs[2] - z)) / (al + bm_plus_cn);
		//x_proj = (fil2->x[0]*al + x_obs*(bm_plus_cn) + fil2->lenvect[1]* fil2->lenvect[0]*(fil1->y[0] - y_obs) + fil2->lenvect[2] * fil2->lenvect[0] * (fil1->z[0] - z_obs)) / (al + bm_plus_cn);
	}
	else
	{
		x_y_z[0] = xyz_obs[0];
	}
	if (grad[1] != 0)
	{
		double al_plus_cn = al + cn;
		x_y_z[1] = (y * bm + xyz_obs[1] * al_plus_cn - grad[0] * grad[1] * (xyz_obs[0] - x) - grad[2] * grad[1] * (xyz_obs[2] - z)) / (al_plus_cn + bm);
	}
	else
	{
		x_y_z[1] = xyz_obs[1];
	}
	if (grad[2] != 0)
	{
		double al_plus_bm = al + bm;
		x_y_z[2] = (z * cn + xyz_obs[2] * (al_plus_bm)-grad[1] * grad[2] * (xyz_obs[1] - y) - grad[0] * grad[2] * (xyz_obs[0] - x)) / (al_plus_bm + cn);
	}
	else
	{
		x_y_z[2] = xyz_obs[2];
	}
}


/*Here the intersection between a line and a plane is calculated*/
double dV_R_temp1(double x_proj, double y_proj, double z_proj, double x_center, double y_center, double z_center, double x_obs, double y_obs, double z_obs, double l[3], double n[3]
	, double delta_center_x, double delta_center_y, double delta_center_z)
{
	double p_plus[3];
	double p_min[3];
	double p0[3];
	double U[3];
	double x_y_z[3];
	double xyz_obs[] = { 0,0,0 };
	xyz_point(x_y_z, n, xyz_obs, x_center, y_center, z_center);
	x_proj = x_proj - x_y_z[0];
	y_proj = y_proj - x_y_z[1];
	z_proj = z_proj - x_y_z[2];
	//U[0] = -U[0];
	//U[1] = -U[1];
	//U[2] = -U[2];
	double n_dot_r = n[0] * (x_center + delta_center_x) + n[1] * (y_center + delta_center_y) + n[2] * (z_center + delta_center_z);
	p_plus[0] = x_center + delta_center_x - n[0] * (n_dot_r);
	p_plus[1] = y_center + delta_center_y - n[1] * (n_dot_r);
	p_plus[2] = z_center + delta_center_z - n[2] * (n_dot_r);


	n_dot_r = n[0] * (x_center - delta_center_x) + n[1] * (y_center - delta_center_y) + n[2] * (z_center - delta_center_z);
	p_min[0] = x_center - delta_center_x - n[0] * (n_dot_r);
	p_min[1] = y_center - delta_center_y - n[1] * (n_dot_r);
	p_min[2] = z_center - delta_center_z - n[2] * (n_dot_r);

	double lplus_ij = (p_plus[0] - x_proj)*l[0] + (p_plus[1] - y_proj)*l[1] + (p_plus[2] - z_proj)*l[2];
	double lmin_ij = (p_min[0] - x_proj)*l[0] + (p_min[1] - y_proj)*l[1] + (p_min[2] - z_proj)*l[2];


	U[0] = l[1] * n[2] - l[2] * n[1];
	U[1] = l[2] * n[0] - l[0] * n[2];
	U[2] = l[0] * n[1] - l[1] * n[0];

	double P0_ij = abs((p_plus[0] - x_proj)*U[0] + (p_plus[1] - y_proj)*U[1] + (p_plus[2] - z_proj)*U[2]);

	if (P0_ij != 0)
	{
		p0[0] = ((p_min[0] - x_proj) - lmin_ij*l[0]) / P0_ij;
		p0[1] = ((p_min[1] - y_proj) - lmin_ij*l[1]) / P0_ij;
		p0[2] = ((p_min[2] - z_proj) - lmin_ij*l[2]) / P0_ij;
	}

	else
	{
		return 0;
		p0[0] = U[0];
		p0[1] = U[1];
		p0[2] = U[2];
	}

	double d_j = abs((x_obs - x_center)*n[0] + (y_obs - y_center)*n[1] + (z_obs - z_center)*n[2]);

	double R0_ij = sqrt(pow(P0_ij, 2) + pow(d_j, 2));
	double p_plus2 = sqrt(pow(P0_ij, 2) + pow(lplus_ij, 2));
	double p_min2 = sqrt(pow(P0_ij, 2) + pow(lmin_ij, 2));
	double Rmin_ij = sqrt(pow(p_min2, 2) + pow(d_j, 2));
	double Rplus_ij = sqrt(pow(p_plus2, 2) + pow(d_j, 2));

	double P0ij_dot_Uij = p0[0] * U[0] + p0[1] * U[1] + p0[2] * U[2];
	if (Rmin_ij == -lmin_ij || Rplus_ij == -lplus_ij) return 0;
	return P0ij_dot_Uij*(d_j*(atan((P0_ij*lplus_ij) / (pow(R0_ij, 2) + d_j*Rplus_ij))
		- atan((P0_ij*lmin_ij) / (pow(R0_ij, 2) + d_j*Rmin_ij)))
		- P0_ij*log((Rplus_ij + lplus_ij) / (Rmin_ij + lmin_ij)));
}


/*Computes Gaussian quadrature integration on the 
outer integral with predetermined number of pointes on width, height and length 
and then analytic integration on the inner integral.
This is typically called to compute self - inductance terms and mutual inductance of near filaments.*/
double aanalitic(FILAMENT *fil1, FILAMENT *fil2, int points, int points2, int points3) {//new points
																						//double a_b = /*fil1->area  * fil1->length*/ 1/ 8;
	int points_length = points;
	int points_width = points2;
	int points_height = points3;






	/*
	if (fil1->length > fil1->width && fil1->length > fil1->height)
	{
	points_length = 5;
	if (fil1->width*4 > fil1->length)
	{
	points_width = 5;
	}
	else
	{
	points_width = 2;
	}
	if (fil1->height * 4 > fil1->length)
	{
	points_height = 5;
	}
	else
	{
	points_height = 2;
	}

	}

	else if (fil1->width > fil1->height)
	{
	points_width = 5;

	if (fil1->length * 4 > fil1->width)
	{
	points_length = 5;
	}
	else
	{
	points_length = 2;
	}
	if (fil1->height * 4 > fil1->width)
	{
	points_height = 5;
	}
	else
	{
	points_height = 2;
	}
	}
	else
	{
	points_height = 5;

	if (fil1->length * 4 > fil1->height)
	{
	points_length = 5;
	}
	else
	{
	points_length = 2;
	}
	if (fil1->width * 4 > fil1->height)
	{
	points_width = 5;
	}
	else
	{
	points_width = 2;
	}

	}
	*/

	double width1[3];
	double height1[3];
	double width2[3];
	double height2[3];
	get_wid(fil1, width1);


	get_height(fil1, width1, height1);

	get_wid(fil2, width2);




	get_height(fil2, width2, height2);



	double *wi = (double *)malloc(points_length * sizeof(double));
	double *wi1 = (double *)malloc(points_width * sizeof(double));
	double *wi2 = (double *)malloc(points_height * sizeof(double));

	double *xi = (double *)malloc(points_length * sizeof(double));
	double *xi1 = (double *)malloc(points_width * sizeof(double));
	double *xi2 = (double *)malloc(points_height * sizeof(double));




	if (points == 1)
	{
		wi[0] = 1.0000000000000000e+000;
	}
	else if (points == 2)
	{
		wi[0] = 4.9999999999999989e-001;
		wi[1] = 4.9999999999999989e-001;
	}
	else if (points == 3)
	{
		wi[0] = 2.7777777777777768e-001;
		wi[1] = 4.4444444444444425e-001;
		wi[2] = 2.7777777777777779e-001;
	}
	else if (points == 4)
	{
		wi[0] = 1.7392742256872706e-001;
		wi[1] = 3.2607257743127316e-001;
		wi[2] = 3.2607257743127299e-001;
		wi[3] = 1.7392742256872701e-001;

	}
	else if (points == 5)
	{
		wi[0] = 1.1846344252809456e-001;
		wi[1] = 2.3931433524968337e-001;
		wi[2] = 2.8444444444444444e-001;
		wi[3] = 2.3931433524968343e-001;
		wi[4] = 1.1846344252809445e-001;
	}
	else if (points == 7)
	{
		wi[0] = 6.4742483084434824e-002;
		wi[1] = 1.3985269574463838e-001;
		wi[2] = 1.9091502525255943e-001;
		wi[3] = 2.0897959183673431e-001;
		wi[4] = 1.9091502525255952e-001;
		wi[5] = 1.3985269574463827e-001;
		wi[6] = 6.4742483084434907e-002;

	}
	else if (points == 8)
	{
		wi[0] = 5.0614268145188081e-002;
		wi[1] = 1.1119051722668732e-001;
		wi[2] = 1.5685332293894375e-001;
		wi[3] = 1.8134189168918130e-001;
		wi[4] = 1.8134189168918088e-001;
		wi[5] = 1.5685332293894350e-001;
		wi[6] = 1.1119051722668724e-001;
		wi[7] = 5.0614268145188233e-002;

	}
	else if (points == 10)
	{
		wi[0] = 3.3335672154344055e-002;
		wi[1] = 7.4725674575290474e-002;
		wi[2] = 1.0954318125799090e-001;
		wi[3] = 1.3463335965499817e-001;
		wi[4] = 1.4776211235737666e-001;
		wi[5] = 1.4776211235737630e-001;
		wi[6] = 1.3463335965499804e-001;
		wi[7] = 1.0954318125799063e-001;
		wi[8] = 7.4725674575290349e-002;
		wi[9] = 3.3335672154344208e-002;

	}
	else if (points == 12)
	{
		wi[0] = 2.3587668193255702e-002;
		wi[1] = 5.3469662997659248e-002;
		wi[2] = 8.0039164271672986e-002;
		wi[3] = 1.0158371336153295e-001;
		wi[4] = 1.1674626826917765e-001;
		wi[5] = 1.2457352290670126e-001;
		wi[6] = 1.2457352290670118e-001;
		wi[7] = 1.1674626826917735e-001;
		wi[8] = 1.0158371336153298e-001;
		wi[9] = 8.0039164271672639e-002;
		wi[10] = 5.3469662997659428e-002;
		wi[11] = 2.3587668193256139e-002;
	}
	else if (points == 13)
	{
		wi[0] = 2.0242002382657838e-002;
		wi[1] = 4.6060749918864448e-002;
		wi[2] = 6.9436755109893736e-002;
		wi[3] = 8.9072990380972800e-002;
		wi[4] = 1.0390802376844431e-001;
		wi[5] = 1.1314159013144860e-001;
		wi[6] = 1.1627577661543718e-001;
		wi[7] = 1.1314159013144845e-001;
		wi[8] = 1.0390802376844471e-001;
		wi[9] = 8.9072990380973188e-002;
		wi[10] = 6.9436755109893292e-002;
		wi[11] = 4.6060749918864434e-002;
		wi[12] = 2.0242002382657821e-002;
	}
	else if (points == 18)
	{
		wi[0] = 1.0808006763241524e-002;
		wi[1] = 2.4857274447485238e-002;
		wi[2] = 3.8212865127444338e-002;
		wi[3] = 5.0471022053143515e-002;
		wi[4] = 6.1277603355739250e-002;
		wi[5] = 7.0321457335325410e-002;
		wi[6] = 7.7342337563132565e-002;
		wi[7] = 8.2138241872916282e-002;
		wi[8] = 8.4571191481571661e-002;
		wi[9] = 8.4571191481571495e-002;
		wi[10] = 8.2138241872916129e-002;
		wi[11] = 7.7342337563132746e-002;
		wi[12] = 7.0321457335324702e-002;
		wi[13] = 6.1277603355739153e-002;
		wi[14] = 5.0471022053143529e-002;
		wi[15] = 3.8212865127444477e-002;
		wi[16] = 2.4857274447485082e-002;
		wi[17] = 1.0808006763241576e-002;
	}
	else if (points == 27)
	{
		wi[0] = 4.8994980256471419e-003;
		wi[1] = 1.1343115798090284e-002;
		wi[2] = 1.7648526878709603e-002;
		wi[3] = 2.3724706260307579e-002;
		wi[4] = 2.9491768429917006e-002;
		wi[5] = 3.4874411883122618e-002;
		wi[6] = 3.9802433886528973e-002;
		wi[7] = 4.4211579271878021e-002;
		wi[8] = 4.8044363685014677e-002;
		wi[9] = 5.1250818908873216e-002;
		wi[10] = 5.3789142894266430e-002;
		wi[11] = 5.5626244178422456e-002;
		wi[12] = 5.6738173054482456e-002;
		wi[13] = 5.7110433689478726e-002;
		wi[14] = 5.6738173054482761e-002;
		wi[15] = 5.5626244178422705e-002;
		wi[16] = 5.3789142894266902e-002;
		wi[17] = 5.1250818908872953e-002;
		wi[18] = 4.8044363685013941e-002;
		wi[19] = 4.4211579271878250e-002;
		wi[20] = 3.9802433886529522e-002;
		wi[21] = 3.4874411883122666e-002;
		wi[22] = 2.9491768429916607e-002;
		wi[23] = 2.3724706260307819e-002;
		wi[24] = 1.7648526878709755e-002;
		wi[25] = 1.1343115798090693e-002;
		wi[26] = 4.8994980256471142e-003;
	}
	else if (points == 31)
	{
		wi[0] = 3.7354157896245144e-003;
		wi[1] = 8.6593103951552244e-003;
		wi[2] = 1.3504509592489613e-002;
		wi[3] = 1.8216136956192775e-002;
		wi[4] = 2.2746853763600677e-002;
		wi[5] = 2.7051541212458532e-002;
		wi[6] = 3.1087393280513912e-002;
		wi[7] = 3.4814291617705204e-002;
		wi[8] = 3.8195193299388357e-002;
		wi[9] = 4.1196495880794611e-002;
		wi[10] = 4.3788370304238822e-002;
		wi[11] = 4.5945056946820530e-002;
		wi[12] = 4.7645121456159768e-002;
		wi[13] = 4.8871667693164013e-002;
		wi[14] = 4.9612505613336480e-002;
		wi[15] = 4.9860272396713541e-002;
		wi[16] = 4.9612505613335800e-002;
		wi[17] = 4.8871667693164256e-002;
		wi[18] = 4.7645121456159671e-002;
		wi[19] = 4.5945056946820877e-002;
		wi[20] = 4.3788370304238822e-002;
		wi[21] = 4.1196495880794437e-002;
		wi[22] = 3.8195193299388988e-002;
		wi[23] = 3.4814291617704482e-002;
		wi[24] = 3.1087393280514110e-002;
		wi[25] = 2.7051541212458823e-002;
		wi[26] = 2.2746853763600399e-002;
		wi[27] = 1.8216136956192744e-002;
		wi[28] = 1.3504509592489819e-002;
		wi[29] = 8.6593103951551568e-003;
		wi[30] = 3.7354157896242767e-003;
	}



	if (points2 == 1)
	{
		wi1[0] = 1.0000000000000000e+000;
	}
	else if (points2 == 2)
	{
		wi1[0] = 4.9999999999999989e-001;
		wi1[1] = 4.9999999999999989e-001;

	}
	else if (points2 == 3)
	{
		wi1[0] = 2.7777777777777768e-001;
		wi1[1] = 4.4444444444444425e-001;
		wi1[2] = 2.7777777777777779e-001;
	}
	else if (points2 == 4)
	{
		wi1[0] = 1.7392742256872706e-001;
		wi1[1] = 3.2607257743127316e-001;
		wi1[2] = 3.2607257743127299e-001;
		wi1[3] = 1.7392742256872701e-001;

	}
	else if (points2 == 5)
	{
		wi1[0] = 1.1846344252809456e-001;
		wi1[1] = 2.3931433524968337e-001;
		wi1[2] = 2.8444444444444444e-001;
		wi1[3] = 2.3931433524968343e-001;
		wi1[4] = 1.1846344252809445e-001;

	}
	else if (points2 == 7)
	{
		wi1[0] = 6.4742483084434824e-002;
		wi1[1] = 1.3985269574463838e-001;
		wi1[2] = 1.9091502525255943e-001;
		wi1[3] = 2.0897959183673431e-001;
		wi1[4] = 1.9091502525255952e-001;
		wi1[5] = 1.3985269574463827e-001;
		wi1[6] = 6.4742483084434907e-002;
	}
	else if (points2 == 8)
	{
		wi1[0] = 5.0614268145188081e-002;
		wi1[1] = 1.1119051722668732e-001;
		wi1[2] = 1.5685332293894375e-001;
		wi1[3] = 1.8134189168918130e-001;
		wi1[4] = 1.8134189168918088e-001;
		wi1[5] = 1.5685332293894350e-001;
		wi1[6] = 1.1119051722668724e-001;
		wi1[7] = 5.0614268145188233e-002;
	}
	else if (points2 == 10)
	{
		wi1[0] = 3.3335672154344055e-002;
		wi1[1] = 7.4725674575290474e-002;
		wi1[2] = 1.0954318125799090e-001;
		wi1[3] = 1.3463335965499817e-001;
		wi1[4] = 1.4776211235737666e-001;
		wi1[5] = 1.4776211235737630e-001;
		wi1[6] = 1.3463335965499804e-001;
		wi1[7] = 1.0954318125799063e-001;
		wi1[8] = 7.4725674575290349e-002;
		wi1[9] = 3.3335672154344208e-002;
	}
	else if (points2 == 12)
	{
		wi1[0] = 2.3587668193255702e-002;
		wi1[1] = 5.3469662997659248e-002;
		wi1[2] = 8.0039164271672986e-002;
		wi1[3] = 1.0158371336153295e-001;
		wi1[4] = 1.1674626826917765e-001;
		wi1[5] = 1.2457352290670126e-001;
		wi1[6] = 1.2457352290670118e-001;
		wi1[7] = 1.1674626826917735e-001;
		wi1[8] = 1.0158371336153298e-001;
		wi1[9] = 8.0039164271672639e-002;
		wi1[10] = 5.3469662997659428e-002;
		wi1[11] = 2.3587668193256139e-002;
	}
	else if (points2 == 13)
	{
		wi1[0] = 2.0242002382657838e-002;
		wi1[1] = 4.6060749918864448e-002;
		wi1[2] = 6.9436755109893736e-002;
		wi1[3] = 8.9072990380972800e-002;
		wi1[4] = 1.0390802376844431e-001;
		wi1[5] = 1.1314159013144860e-001;
		wi1[6] = 1.1627577661543718e-001;
		wi1[7] = 1.1314159013144845e-001;
		wi1[8] = 1.0390802376844471e-001;
		wi1[9] = 8.9072990380973188e-002;
		wi1[10] = 6.9436755109893292e-002;
		wi1[11] = 4.6060749918864434e-002;
		wi1[12] = 2.0242002382657821e-002;
	}
	else if (points2 == 18)
	{
		wi1[0] = 1.0808006763241524e-002;
		wi1[1] = 2.4857274447485238e-002;
		wi1[2] = 3.8212865127444338e-002;
		wi1[3] = 5.0471022053143515e-002;
		wi1[4] = 6.1277603355739250e-002;
		wi1[5] = 7.0321457335325410e-002;
		wi1[6] = 7.7342337563132565e-002;
		wi1[7] = 8.2138241872916282e-002;
		wi1[8] = 8.4571191481571661e-002;
		wi1[9] = 8.4571191481571495e-002;
		wi1[10] = 8.2138241872916129e-002;
		wi1[11] = 7.7342337563132746e-002;
		wi1[12] = 7.0321457335324702e-002;
		wi1[13] = 6.1277603355739153e-002;
		wi1[14] = 5.0471022053143529e-002;
		wi1[15] = 3.8212865127444477e-002;
		wi1[16] = 2.4857274447485082e-002;
		wi1[17] = 1.0808006763241576e-002;
	}
	else if (points2 == 27)
	{
		wi1[0] = 4.8994980256471419e-003;
		wi1[1] = 1.1343115798090284e-002;
		wi1[2] = 1.7648526878709603e-002;
		wi1[3] = 2.3724706260307579e-002;
		wi1[4] = 2.9491768429917006e-002;
		wi1[5] = 3.4874411883122618e-002;
		wi1[6] = 3.9802433886528973e-002;
		wi1[7] = 4.4211579271878021e-002;
		wi1[8] = 4.8044363685014677e-002;
		wi1[9] = 5.1250818908873216e-002;
		wi1[10] = 5.3789142894266430e-002;
		wi1[11] = 5.5626244178422456e-002;
		wi1[12] = 5.6738173054482456e-002;
		wi1[13] = 5.7110433689478726e-002;
		wi1[14] = 5.6738173054482761e-002;
		wi1[15] = 5.5626244178422705e-002;
		wi1[16] = 5.3789142894266902e-002;
		wi1[17] = 5.1250818908872953e-002;
		wi1[18] = 4.8044363685013941e-002;
		wi1[19] = 4.4211579271878250e-002;
		wi1[20] = 3.9802433886529522e-002;
		wi1[21] = 3.4874411883122666e-002;
		wi1[22] = 2.9491768429916607e-002;
		wi1[23] = 2.3724706260307819e-002;
		wi1[24] = 1.7648526878709755e-002;
		wi1[25] = 1.1343115798090693e-002;
		wi1[26] = 4.8994980256471142e-003;
	}
	else if (points2 == 31)
	{
		wi1[0] = 3.7354157896245144e-003;
		wi1[1] = 8.6593103951552244e-003;
		wi1[2] = 1.3504509592489613e-002;
		wi1[3] = 1.8216136956192775e-002;
		wi1[4] = 2.2746853763600677e-002;
		wi1[5] = 2.7051541212458532e-002;
		wi1[6] = 3.1087393280513912e-002;
		wi1[7] = 3.4814291617705204e-002;
		wi1[8] = 3.8195193299388357e-002;
		wi1[9] = 4.1196495880794611e-002;
		wi1[10] = 4.3788370304238822e-002;
		wi1[11] = 4.5945056946820530e-002;
		wi1[12] = 4.7645121456159768e-002;
		wi1[13] = 4.8871667693164013e-002;
		wi1[14] = 4.9612505613336480e-002;
		wi1[15] = 4.9860272396713541e-002;
		wi1[16] = 4.9612505613335800e-002;
		wi1[17] = 4.8871667693164256e-002;
		wi1[18] = 4.7645121456159671e-002;
		wi1[19] = 4.5945056946820877e-002;
		wi1[20] = 4.3788370304238822e-002;
		wi1[21] = 4.1196495880794437e-002;
		wi1[22] = 3.8195193299388988e-002;
		wi1[23] = 3.4814291617704482e-002;
		wi1[24] = 3.1087393280514110e-002;
		wi1[25] = 2.7051541212458823e-002;
		wi1[26] = 2.2746853763600399e-002;
		wi1[27] = 1.8216136956192744e-002;
		wi1[28] = 1.3504509592489819e-002;
		wi1[29] = 8.6593103951551568e-003;
		wi1[30] = 3.7354157896242767e-003;
	}



	if (points3 == 1)
	{
		wi2[0] = 1.0000000000000000e+000;
	}
	else if (points3 == 2)
	{
		wi2[0] = 4.9999999999999989e-001;
		wi2[1] = 4.9999999999999989e-001;
	}
	else if (points3 == 3)
	{
		wi2[0] = 2.7777777777777768e-001;
		wi2[1] = 4.4444444444444425e-001;
		wi2[2] = 2.7777777777777779e-001;
	}
	else if (points3 == 4)
	{
		wi2[0] = 1.7392742256872706e-001;
		wi2[1] = 3.2607257743127316e-001;
		wi2[2] = 3.2607257743127299e-001;
		wi2[3] = 1.7392742256872701e-001;
	}
	else if (points3 == 5)
	{
		wi2[0] = 1.1846344252809456e-001;
		wi2[1] = 2.3931433524968337e-001;
		wi2[2] = 2.8444444444444444e-001;
		wi2[3] = 2.3931433524968343e-001;
		wi2[4] = 1.1846344252809445e-001;
	}
	else if (points3 == 7)
	{
		wi2[0] = 6.4742483084434824e-002;
		wi2[1] = 1.3985269574463838e-001;
		wi2[2] = 1.9091502525255943e-001;
		wi2[3] = 2.0897959183673431e-001;
		wi2[4] = 1.9091502525255952e-001;
		wi2[5] = 1.3985269574463827e-001;
		wi2[6] = 6.4742483084434907e-002;
	}
	else if (points3 == 8)
	{
		wi2[0] = 5.0614268145188081e-002;
		wi2[1] = 1.1119051722668732e-001;
		wi2[2] = 1.5685332293894375e-001;
		wi2[3] = 1.8134189168918130e-001;
		wi2[4] = 1.8134189168918088e-001;
		wi2[5] = 1.5685332293894350e-001;
		wi2[6] = 1.1119051722668724e-001;
		wi2[7] = 5.0614268145188233e-002;
	}
	else if (points3 == 10)
	{
		wi2[0] = 3.3335672154344055e-002;
		wi2[1] = 7.4725674575290474e-002;
		wi2[2] = 1.0954318125799090e-001;
		wi2[3] = 1.3463335965499817e-001;
		wi2[4] = 1.4776211235737666e-001;
		wi2[5] = 1.4776211235737630e-001;
		wi2[6] = 1.3463335965499804e-001;
		wi2[7] = 1.0954318125799063e-001;
		wi2[8] = 7.4725674575290349e-002;
		wi2[9] = 3.3335672154344208e-002;
	}
	else if (points3 == 12)
	{
		wi2[0] = 2.3587668193255702e-002;
		wi2[1] = 5.3469662997659248e-002;
		wi2[2] = 8.0039164271672986e-002;
		wi2[3] = 1.0158371336153295e-001;
		wi2[4] = 1.1674626826917765e-001;
		wi2[5] = 1.2457352290670126e-001;
		wi2[6] = 1.2457352290670118e-001;
		wi2[7] = 1.1674626826917735e-001;
		wi2[8] = 1.0158371336153298e-001;
		wi2[9] = 8.0039164271672639e-002;
		wi2[10] = 5.3469662997659428e-002;
		wi2[11] = 2.3587668193256139e-002;
	}
	else if (points3 == 13)
	{
		wi2[0] = 2.0242002382657838e-002;
		wi2[1] = 4.6060749918864448e-002;
		wi2[2] = 6.9436755109893736e-002;
		wi2[3] = 8.9072990380972800e-002;
		wi2[4] = 1.0390802376844431e-001;
		wi2[5] = 1.1314159013144860e-001;
		wi2[6] = 1.1627577661543718e-001;
		wi2[7] = 1.1314159013144845e-001;
		wi2[8] = 1.0390802376844471e-001;
		wi2[9] = 8.9072990380973188e-002;
		wi2[10] = 6.9436755109893292e-002;
		wi2[11] = 4.6060749918864434e-002;
		wi2[12] = 2.0242002382657821e-002;
	}
	else if (points3 == 18)
	{
		wi2[0] = 1.0808006763241524e-002;
		wi2[1] = 2.4857274447485238e-002;
		wi2[2] = 3.8212865127444338e-002;
		wi2[3] = 5.0471022053143515e-002;
		wi2[4] = 6.1277603355739250e-002;
		wi2[5] = 7.0321457335325410e-002;
		wi2[6] = 7.7342337563132565e-002;
		wi2[7] = 8.2138241872916282e-002;
		wi2[8] = 8.4571191481571661e-002;
		wi2[9] = 8.4571191481571495e-002;
		wi2[10] = 8.2138241872916129e-002;
		wi2[11] = 7.7342337563132746e-002;
		wi2[12] = 7.0321457335324702e-002;
		wi2[13] = 6.1277603355739153e-002;
		wi2[14] = 5.0471022053143529e-002;
		wi2[15] = 3.8212865127444477e-002;
		wi2[16] = 2.4857274447485082e-002;
		wi2[17] = 1.0808006763241576e-002;
	}
	else if (points3 == 27)
	{
		wi2[0] = 4.8994980256471419e-003;
		wi2[1] = 1.1343115798090284e-002;
		wi2[2] = 1.7648526878709603e-002;
		wi2[3] = 2.3724706260307579e-002;
		wi2[4] = 2.9491768429917006e-002;
		wi2[5] = 3.4874411883122618e-002;
		wi2[6] = 3.9802433886528973e-002;
		wi2[7] = 4.4211579271878021e-002;
		wi2[8] = 4.8044363685014677e-002;
		wi2[9] = 5.1250818908873216e-002;
		wi2[10] = 5.3789142894266430e-002;
		wi2[11] = 5.5626244178422456e-002;
		wi2[12] = 5.6738173054482456e-002;
		wi2[13] = 5.7110433689478726e-002;
		wi2[14] = 5.6738173054482761e-002;
		wi2[15] = 5.5626244178422705e-002;
		wi2[16] = 5.3789142894266902e-002;
		wi2[17] = 5.1250818908872953e-002;
		wi2[18] = 4.8044363685013941e-002;
		wi2[19] = 4.4211579271878250e-002;
		wi2[20] = 3.9802433886529522e-002;
		wi2[21] = 3.4874411883122666e-002;
		wi2[22] = 2.9491768429916607e-002;
		wi2[23] = 2.3724706260307819e-002;
		wi2[24] = 1.7648526878709755e-002;
		wi2[25] = 1.1343115798090693e-002;
		wi2[26] = 4.8994980256471142e-003;
	}
	else if (points3 == 31)
	{
		wi2[0] = 3.7354157896245144e-003;
		wi2[1] = 8.6593103951552244e-003;
		wi2[2] = 1.3504509592489613e-002;
		wi2[3] = 1.8216136956192775e-002;
		wi2[4] = 2.2746853763600677e-002;
		wi2[5] = 2.7051541212458532e-002;
		wi2[6] = 3.1087393280513912e-002;
		wi2[7] = 3.4814291617705204e-002;
		wi2[8] = 3.8195193299388357e-002;
		wi2[9] = 4.1196495880794611e-002;
		wi2[10] = 4.3788370304238822e-002;
		wi2[11] = 4.5945056946820530e-002;
		wi2[12] = 4.7645121456159768e-002;
		wi2[13] = 4.8871667693164013e-002;
		wi2[14] = 4.9612505613336480e-002;
		wi2[15] = 4.9860272396713541e-002;
		wi2[16] = 4.9612505613335800e-002;
		wi2[17] = 4.8871667693164256e-002;
		wi2[18] = 4.7645121456159671e-002;
		wi2[19] = 4.5945056946820877e-002;
		wi2[20] = 4.3788370304238822e-002;
		wi2[21] = 4.1196495880794437e-002;
		wi2[22] = 3.8195193299388988e-002;
		wi2[23] = 3.4814291617704482e-002;
		wi2[24] = 3.1087393280514110e-002;
		wi2[25] = 2.7051541212458823e-002;
		wi2[26] = 2.2746853763600399e-002;
		wi2[27] = 1.8216136956192744e-002;
		wi2[28] = 1.3504509592489819e-002;
		wi2[29] = 8.6593103951551568e-003;
		wi2[30] = 3.7354157896242767e-003;
	}



	for (int i = 0; i < points; i++)
	{
		wi[i] *= 2;
	}
	for (int i = 0; i < points2; i++)
	{
		wi1[i] *= 2;
	}
	for (int i = 0; i < points3; i++)
	{
		wi2[i] *= 2;
	}


	if (points == 1)
	{
		xi[0] = 5.0000000000000000e-001;
	}
	else if (points == 2)
	{

		xi[0] = 2.1132486540518708e-001;
		xi[1] = 7.8867513459481287e-001;

	}
	else if (points == 3)
	{

		xi[0] = 1.1270166537925824e-001;
		xi[1] = 5.0000000000000000e-001;
		xi[2] = 8.8729833462074170e-001;
	}
	else if (points == 4)
	{

		xi[0] = 6.9431844202973714e-002;
		xi[1] = 3.3000947820757182e-001;
		xi[2] = 6.6999052179242813e-001;
		xi[3] = 9.3056815579702634e-001;
	}
	else if (points == 5)
	{
		xi[0] = 4.6910077030668074e-002;
		xi[1] = 2.3076534494715845e-001;
		xi[2] = 5.0000000000000000e-001;
		xi[3] = 7.6923465505284150e-001;
		xi[4] = 9.5308992296933193e-001;
	}
	else if (points == 7)
	{

		xi[0] = 2.5446043828620812e-002;
		xi[1] = 1.2923440720030277e-001;
		xi[2] = 2.9707742431130146e-001;
		xi[3] = 4.9999999999999994e-001;
		xi[4] = 7.0292257568869854e-001;
		xi[5] = 8.7076559279969734e-001;
		xi[6] = 9.7455395617137919e-001;
	}
	else if (points == 8)
	{

		xi[0] = 1.9855071751231967e-002;
		xi[1] = 1.0166676129318664e-001;
		xi[2] = 2.3723379504183539e-001;
		xi[3] = 4.0828267875217511e-001;
		xi[4] = 5.9171732124782506e-001;
		xi[5] = 7.6276620495816438e-001;
		xi[6] = 8.9833323870681348e-001;
		xi[7] = 9.8014492824876820e-001;
	}
	else if (points == 10)
	{

		xi[0] = 1.3046735741414073e-002;
		xi[1] = 6.7468316655507898e-002;
		xi[2] = 1.6029521585048767e-001;
		xi[3] = 2.8330230293537639e-001;
		xi[4] = 4.2556283050918453e-001;
		xi[5] = 5.7443716949081558e-001;
		xi[6] = 7.1669769706462361e-001;
		xi[7] = 8.3970478414951222e-001;
		xi[8] = 9.3253168334449221e-001;
		xi[9] = 9.8695326425858587e-001;
	}
	else if (points == 12)
	{
		xi[0] = 9.2196828766404337e-003;
		xi[1] = 4.7941371814762324e-002;
		xi[2] = 1.1504866290284754e-001;
		xi[3] = 2.0634102285669120e-001;
		xi[4] = 3.1608425050090982e-001;
		xi[5] = 4.3738329574426549e-001;
		xi[6] = 5.6261670425573440e-001;
		xi[7] = 6.8391574949909018e-001;
		xi[8] = 7.9365897714330869e-001;
		xi[9] = 8.8495133709715246e-001;
		xi[10] = 9.5205862818523768e-001;
		xi[11] = 9.9078031712335957e-001;
	}
	else if (points == 13)
	{
		xi[0] = 7.9084726407059325e-003;
		xi[1] = 4.1200800388511094e-002;
		xi[2] = 9.9210954633345005e-002;
		xi[3] = 1.7882533027982983e-001;
		xi[4] = 2.7575362448177665e-001;
		xi[5] = 3.8477084202243278e-001;
		xi[6] = 4.9999999999999989e-001;
		xi[7] = 6.1522915797756739e-001;
		xi[8] = 7.2424637551822357e-001;
		xi[9] = 8.2117466972017006e-001;
		xi[10] = 9.0078904536665494e-001;
		xi[11] = 9.5879919961148885e-001;
		xi[12] = 9.9209152735929407e-001;
	}
	else if (points == 18)
	{
		xi[0] = 4.2174157895344955e-003;
		xi[1] = 2.2088025214301199e-002;
		xi[2] = 5.3698766751222149e-002;
		xi[3] = 9.8147520513738540e-002;
		xi[4] = 1.5415647846982339e-001;
		xi[5] = 2.2011458446302623e-001;
		xi[6] = 2.9412441926857869e-001;
		xi[7] = 3.7405688715424718e-001;
		xi[8] = 4.5761249347913252e-001;
		xi[9] = 5.4238750652086776e-001;
		xi[10] = 6.2594311284575288e-001;
		xi[11] = 7.0587558073142143e-001;
		xi[12] = 7.7988541553697388e-001;
		xi[13] = 8.4584352153017661e-001;
		xi[14] = 9.0185247948626157e-001;
		xi[15] = 9.4630123324877791e-001;
		xi[16] = 9.7791197478569891e-001;
		xi[17] = 9.9578258421046550e-001;
	}
	else if (points == 27)
	{
		xi[0] = 1.9103685555057481e-003;
		xi[1] = 1.0038262019249400e-002;
		xi[2] = 2.4549721092647525e-002;
		xi[3] = 4.5258839661254380e-002;
		xi[4] = 7.1896045990852642e-002;
		xi[5] = 1.0411418046474580e-001;
		xi[6] = 1.4149326313028820e-001;
		xi[7] = 1.8354601402675236e-001;
		xi[8] = 2.2972421771027135e-001;
		xi[9] = 2.7942587412498643e-001;
		xi[10] = 3.3200304818074544e-001;
		xi[11] = 3.8677031728023159e-001;
		xi[12] = 4.4301370719523503e-001;
		xi[13] = 5.0000000000000000e-001;
		xi[14] = 5.5698629280476519e-001;
		xi[15] = 6.1322968271976863e-001;
		xi[16] = 6.6799695181925434e-001;
		xi[17] = 7.2057412587501357e-001;
		xi[18] = 7.7027578228972848e-001;
		xi[19] = 8.1645398597324770e-001;
		xi[20] = 8.5850673686971191e-001;
		xi[21] = 8.9588581953525415e-001;
		xi[22] = 9.2810395400914714e-001;
		xi[23] = 9.5474116033874556e-001;
		xi[24] = 9.7545027890735247e-001;
		xi[25] = 9.8996173798075060e-001;
		xi[26] = 9.9808963144449425e-001;
	}
	else if (points == 31)
	{
		xi[0] = 1.4562590902615358e-003;
		xi[1] = 7.6570451674239925e-003;
		xi[2] = 1.8748037453524935e-002;
		xi[3] = 3.4621501051675774e-002;
		xi[4] = 5.5119985025864404e-002;
		xi[5] = 8.0039839926866363e-002;
		xi[6] = 1.0913342579168744e-001;
		xi[7] = 1.4211160770657344e-001;
		xi[8] = 1.7864663853786994e-001;
		xi[9] = 2.1837541929642545e-001;
		xi[10] = 2.6090310897754887e-001;
		xi[11] = 3.0580704919588364e-001;
		xi[12] = 3.5264096500914910e-001;
		xi[13] = 4.0093940033221465e-001;
		xi[14] = 4.5022234392382915e-001;
		xi[15] = 5.0000000000000000e-001;
		xi[16] = 5.4977765607617068e-001;
		xi[17] = 5.9906059966778546e-001;
		xi[18] = 6.4735903499085090e-001;
		xi[19] = 6.9419295080411647e-001;
		xi[20] = 7.3909689102245146e-001;
		xi[21] = 7.8162458070357455e-001;
		xi[22] = 8.2135336146213023e-001;
		xi[23] = 8.5788839229342662e-001;
		xi[24] = 8.9086657420831239e-001;
		xi[25] = 9.1996016007313386e-001;
		xi[26] = 9.4488001497413532e-001;
		xi[27] = 9.6537849894832406e-001;
		xi[28] = 9.8125196254647484e-001;
		xi[29] = 9.9234295483257617e-001;
		xi[30] = 9.9854374090973852e-001;
	}




	if (points2 == 1)
	{
		xi1[0] = 5.0000000000000000e-001;
	}
	else if (points2 == 2)
	{
		xi1[0] = 2.1132486540518708e-001;
		xi1[1] = 7.8867513459481287e-001;
	}
	else if (points2 == 3)
	{
		xi1[0] = 1.1270166537925824e-001;
		xi1[1] = 5.0000000000000000e-001;
		xi1[2] = 8.8729833462074170e-001;
	}
	else if (points2 == 4)
	{
		xi1[0] = 6.9431844202973714e-002;
		xi1[1] = 3.3000947820757182e-001;
		xi1[2] = 6.6999052179242813e-001;
		xi1[3] = 9.3056815579702634e-001;
	}
	else if (points2 == 5)
	{
		xi1[0] = 4.6910077030668074e-002;
		xi1[1] = 2.3076534494715845e-001;
		xi1[2] = 5.0000000000000000e-001;
		xi1[3] = 7.6923465505284150e-001;
		xi1[4] = 9.5308992296933193e-001;
	}
	else if (points2 == 7)
	{
		xi1[0] = 2.5446043828620812e-002;
		xi1[1] = 1.2923440720030277e-001;
		xi1[2] = 2.9707742431130146e-001;
		xi1[3] = 4.9999999999999994e-001;
		xi1[4] = 7.0292257568869854e-001;
		xi1[5] = 8.7076559279969734e-001;
		xi1[6] = 9.7455395617137919e-001;
	}
	else if (points2 == 8)
	{
		xi1[0] = 1.9855071751231967e-002;
		xi1[1] = 1.0166676129318664e-001;
		xi1[2] = 2.3723379504183539e-001;
		xi1[3] = 4.0828267875217511e-001;
		xi1[4] = 5.9171732124782506e-001;
		xi1[5] = 7.6276620495816438e-001;
		xi1[6] = 8.9833323870681348e-001;
		xi1[7] = 9.8014492824876820e-001;
	}
	else if (points2 == 10)
	{
		xi1[0] = 1.3046735741414073e-002;
		xi1[1] = 6.7468316655507898e-002;
		xi1[2] = 1.6029521585048767e-001;
		xi1[3] = 2.8330230293537639e-001;
		xi1[4] = 4.2556283050918453e-001;
		xi1[5] = 5.7443716949081558e-001;
		xi1[6] = 7.1669769706462361e-001;
		xi1[7] = 8.3970478414951222e-001;
		xi1[8] = 9.3253168334449221e-001;
		xi1[9] = 9.8695326425858587e-001;
	}
	else if (points2 == 12)
	{
		xi1[0] = 9.2196828766404337e-003;
		xi1[1] = 4.7941371814762324e-002;
		xi1[2] = 1.1504866290284754e-001;
		xi1[3] = 2.0634102285669120e-001;
		xi1[4] = 3.1608425050090982e-001;
		xi1[5] = 4.3738329574426549e-001;
		xi1[6] = 5.6261670425573440e-001;
		xi1[7] = 6.8391574949909018e-001;
		xi1[8] = 7.9365897714330869e-001;
		xi1[9] = 8.8495133709715246e-001;
		xi1[10] = 9.5205862818523768e-001;
		xi1[11] = 9.9078031712335957e-001;

	}
	else if (points2 == 13)
	{
		xi1[0] = 7.9084726407059325e-003;
		xi1[1] = 4.1200800388511094e-002;
		xi1[2] = 9.9210954633345005e-002;
		xi1[3] = 1.7882533027982983e-001;
		xi1[4] = 2.7575362448177665e-001;
		xi1[5] = 3.8477084202243278e-001;
		xi1[6] = 4.9999999999999989e-001;
		xi1[7] = 6.1522915797756739e-001;
		xi1[8] = 7.2424637551822357e-001;
		xi1[9] = 8.2117466972017006e-001;
		xi1[10] = 9.0078904536665494e-001;
		xi1[11] = 9.5879919961148885e-001;
		xi1[12] = 9.9209152735929407e-001;
	}
	else if (points2 == 18)
	{
		xi1[0] = 4.2174157895344955e-003;
		xi1[1] = 2.2088025214301199e-002;
		xi1[2] = 5.3698766751222149e-002;
		xi1[3] = 9.8147520513738540e-002;
		xi1[4] = 1.5415647846982339e-001;
		xi1[5] = 2.2011458446302623e-001;
		xi1[6] = 2.9412441926857869e-001;
		xi1[7] = 3.7405688715424718e-001;
		xi1[8] = 4.5761249347913252e-001;
		xi1[9] = 5.4238750652086776e-001;
		xi1[10] = 6.2594311284575288e-001;
		xi1[11] = 7.0587558073142143e-001;
		xi1[12] = 7.7988541553697388e-001;
		xi1[13] = 8.4584352153017661e-001;
		xi1[14] = 9.0185247948626157e-001;
		xi1[15] = 9.4630123324877791e-001;
		xi1[16] = 9.7791197478569891e-001;
		xi1[17] = 9.9578258421046550e-001;
	}
	else if (points2 == 27)
	{
		xi1[0] = 1.9103685555057481e-003;
		xi1[1] = 1.0038262019249400e-002;
		xi1[2] = 2.4549721092647525e-002;
		xi1[3] = 4.5258839661254380e-002;
		xi1[4] = 7.1896045990852642e-002;
		xi1[5] = 1.0411418046474580e-001;
		xi1[6] = 1.4149326313028820e-001;
		xi1[7] = 1.8354601402675236e-001;
		xi1[8] = 2.2972421771027135e-001;
		xi1[9] = 2.7942587412498643e-001;
		xi1[10] = 3.3200304818074544e-001;
		xi1[11] = 3.8677031728023159e-001;
		xi1[12] = 4.4301370719523503e-001;
		xi1[13] = 5.0000000000000000e-001;
		xi1[14] = 5.5698629280476519e-001;
		xi1[15] = 6.1322968271976863e-001;
		xi1[16] = 6.6799695181925434e-001;
		xi1[17] = 7.2057412587501357e-001;
		xi1[18] = 7.7027578228972848e-001;
		xi1[19] = 8.1645398597324770e-001;
		xi1[20] = 8.5850673686971191e-001;
		xi1[21] = 8.9588581953525415e-001;
		xi1[22] = 9.2810395400914714e-001;
		xi1[23] = 9.5474116033874556e-001;
		xi1[24] = 9.7545027890735247e-001;
		xi1[25] = 9.8996173798075060e-001;
		xi1[26] = 9.9808963144449425e-001;
	}
	else if (points2 == 31)
	{
		xi1[0] = 1.4562590902615358e-003;
		xi1[1] = 7.6570451674239925e-003;
		xi1[2] = 1.8748037453524935e-002;
		xi1[3] = 3.4621501051675774e-002;
		xi1[4] = 5.5119985025864404e-002;
		xi1[5] = 8.0039839926866363e-002;
		xi1[6] = 1.0913342579168744e-001;
		xi1[7] = 1.4211160770657344e-001;
		xi1[8] = 1.7864663853786994e-001;
		xi1[9] = 2.1837541929642545e-001;
		xi1[10] = 2.6090310897754887e-001;
		xi1[11] = 3.0580704919588364e-001;
		xi1[12] = 3.5264096500914910e-001;
		xi1[13] = 4.0093940033221465e-001;
		xi1[14] = 4.5022234392382915e-001;
		xi1[15] = 5.0000000000000000e-001;
		xi1[16] = 5.4977765607617068e-001;
		xi1[17] = 5.9906059966778546e-001;
		xi1[18] = 6.4735903499085090e-001;
		xi1[19] = 6.9419295080411647e-001;
		xi1[20] = 7.3909689102245146e-001;
		xi1[21] = 7.8162458070357455e-001;
		xi1[22] = 8.2135336146213023e-001;
		xi1[23] = 8.5788839229342662e-001;
		xi1[24] = 8.9086657420831239e-001;
		xi1[25] = 9.1996016007313386e-001;
		xi1[26] = 9.4488001497413532e-001;
		xi1[27] = 9.6537849894832406e-001;
		xi1[28] = 9.8125196254647484e-001;
		xi1[29] = 9.9234295483257617e-001;
		xi1[30] = 9.9854374090973852e-001;
	}



	if (points3 == 1)
	{
		xi2[0] = 5.0000000000000000e-001;
	}
	else if (points3 == 2)
	{
		xi2[0] = 2.1132486540518708e-001;
		xi2[1] = 7.8867513459481287e-001;
	}
	else if (points3 == 3)
	{
		xi2[0] = 1.1270166537925824e-001;
		xi2[1] = 5.0000000000000000e-001;
		xi2[2] = 8.8729833462074170e-001;
	}
	else if (points3 == 4)
	{
		xi2[0] = 6.9431844202973714e-002;
		xi2[1] = 3.3000947820757182e-001;
		xi2[2] = 6.6999052179242813e-001;
		xi2[3] = 9.3056815579702634e-001;
	}
	else if (points3 == 5)
	{
		xi2[0] = 4.6910077030668074e-002;
		xi2[1] = 2.3076534494715845e-001;
		xi2[2] = 5.0000000000000000e-001;
		xi2[3] = 7.6923465505284150e-001;
		xi2[4] = 9.5308992296933193e-001;
	}
	else if (points3 == 7)
	{
		xi2[0] = 2.5446043828620812e-002;
		xi2[1] = 1.2923440720030277e-001;
		xi2[2] = 2.9707742431130146e-001;
		xi2[3] = 4.9999999999999994e-001;
		xi2[4] = 7.0292257568869854e-001;
		xi2[5] = 8.7076559279969734e-001;
		xi2[6] = 9.7455395617137919e-001;
	}
	else if (points3 == 8)
	{
		xi2[0] = 1.9855071751231967e-002;
		xi2[1] = 1.0166676129318664e-001;
		xi2[2] = 2.3723379504183539e-001;
		xi2[3] = 4.0828267875217511e-001;
		xi2[4] = 5.9171732124782506e-001;
		xi2[5] = 7.6276620495816438e-001;
		xi2[6] = 8.9833323870681348e-001;
		xi2[7] = 9.8014492824876820e-001;
	}
	else if (points3 == 10)
	{
		xi2[0] = 1.3046735741414073e-002;
		xi2[1] = 6.7468316655507898e-002;
		xi2[2] = 1.6029521585048767e-001;
		xi2[3] = 2.8330230293537639e-001;
		xi2[4] = 4.2556283050918453e-001;
		xi2[5] = 5.7443716949081558e-001;
		xi2[6] = 7.1669769706462361e-001;
		xi2[7] = 8.3970478414951222e-001;
		xi2[8] = 9.3253168334449221e-001;
		xi2[9] = 9.8695326425858587e-001;
	}
	else if (points3 == 12)
	{
		xi2[0] = 9.2196828766404337e-003;
		xi2[1] = 4.7941371814762324e-002;
		xi2[2] = 1.1504866290284754e-001;
		xi2[3] = 2.0634102285669120e-001;
		xi2[4] = 3.1608425050090982e-001;
		xi2[5] = 4.3738329574426549e-001;
		xi2[6] = 5.6261670425573440e-001;
		xi2[7] = 6.8391574949909018e-001;
		xi2[8] = 7.9365897714330869e-001;
		xi2[9] = 8.8495133709715246e-001;
		xi2[10] = 9.5205862818523768e-001;
		xi2[11] = 9.9078031712335957e-001;
	}
	else if (points3 == 13)
	{
		xi2[0] = 7.9084726407059325e-003;
		xi2[1] = 4.1200800388511094e-002;
		xi2[2] = 9.9210954633345005e-002;
		xi2[3] = 1.7882533027982983e-001;
		xi2[4] = 2.7575362448177665e-001;
		xi2[5] = 3.8477084202243278e-001;
		xi2[6] = 4.9999999999999989e-001;
		xi2[7] = 6.1522915797756739e-001;
		xi2[8] = 7.2424637551822357e-001;
		xi2[9] = 8.2117466972017006e-001;
		xi2[10] = 9.0078904536665494e-001;
		xi2[11] = 9.5879919961148885e-001;
		xi2[12] = 9.9209152735929407e-001;
	}
	else if (points3 == 18)
	{
		xi2[0] = 4.2174157895344955e-003;
		xi2[1] = 2.2088025214301199e-002;
		xi2[2] = 5.3698766751222149e-002;
		xi2[3] = 9.8147520513738540e-002;
		xi2[4] = 1.5415647846982339e-001;
		xi2[5] = 2.2011458446302623e-001;
		xi2[6] = 2.9412441926857869e-001;
		xi2[7] = 3.7405688715424718e-001;
		xi2[8] = 4.5761249347913252e-001;
		xi2[9] = 5.4238750652086776e-001;
		xi2[10] = 6.2594311284575288e-001;
		xi2[11] = 7.0587558073142143e-001;
		xi2[12] = 7.7988541553697388e-001;
		xi2[13] = 8.4584352153017661e-001;
		xi2[14] = 9.0185247948626157e-001;
		xi2[15] = 9.4630123324877791e-001;
		xi2[16] = 9.7791197478569891e-001;
		xi2[17] = 9.9578258421046550e-001;
	}
	else if (points3 == 27)
	{
		xi2[0] = 1.9103685555057481e-003;
		xi2[1] = 1.0038262019249400e-002;
		xi2[2] = 2.4549721092647525e-002;
		xi2[3] = 4.5258839661254380e-002;
		xi2[4] = 7.1896045990852642e-002;
		xi2[5] = 1.0411418046474580e-001;
		xi2[6] = 1.4149326313028820e-001;
		xi2[7] = 1.8354601402675236e-001;
		xi2[8] = 2.2972421771027135e-001;
		xi2[9] = 2.7942587412498643e-001;
		xi2[10] = 3.3200304818074544e-001;
		xi2[11] = 3.8677031728023159e-001;
		xi2[12] = 4.4301370719523503e-001;
		xi2[13] = 5.0000000000000000e-001;
		xi2[14] = 5.5698629280476519e-001;
		xi2[15] = 6.1322968271976863e-001;
		xi2[16] = 6.6799695181925434e-001;
		xi2[17] = 7.2057412587501357e-001;
		xi2[18] = 7.7027578228972848e-001;
		xi2[19] = 8.1645398597324770e-001;
		xi2[20] = 8.5850673686971191e-001;
		xi2[21] = 8.9588581953525415e-001;
		xi2[22] = 9.2810395400914714e-001;
		xi2[23] = 9.5474116033874556e-001;
		xi2[24] = 9.7545027890735247e-001;
		xi2[25] = 9.8996173798075060e-001;
		xi2[26] = 9.9808963144449425e-001;
	}
	else if (points3 == 31)
	{
		xi2[0] = 1.4562590902615358e-003;
		xi2[1] = 7.6570451674239925e-003;
		xi2[2] = 1.8748037453524935e-002;
		xi2[3] = 3.4621501051675774e-002;
		xi2[4] = 5.5119985025864404e-002;
		xi2[5] = 8.0039839926866363e-002;
		xi2[6] = 1.0913342579168744e-001;
		xi2[7] = 1.4211160770657344e-001;
		xi2[8] = 1.7864663853786994e-001;
		xi2[9] = 2.1837541929642545e-001;
		xi2[10] = 2.6090310897754887e-001;
		xi2[11] = 3.0580704919588364e-001;
		xi2[12] = 3.5264096500914910e-001;
		xi2[13] = 4.0093940033221465e-001;
		xi2[14] = 4.5022234392382915e-001;
		xi2[15] = 5.0000000000000000e-001;
		xi2[16] = 5.4977765607617068e-001;
		xi2[17] = 5.9906059966778546e-001;
		xi2[18] = 6.4735903499085090e-001;
		xi2[19] = 6.9419295080411647e-001;
		xi2[20] = 7.3909689102245146e-001;
		xi2[21] = 7.8162458070357455e-001;
		xi2[22] = 8.2135336146213023e-001;
		xi2[23] = 8.5788839229342662e-001;
		xi2[24] = 8.9086657420831239e-001;
		xi2[25] = 9.1996016007313386e-001;
		xi2[26] = 9.4488001497413532e-001;
		xi2[27] = 9.6537849894832406e-001;
		xi2[28] = 9.8125196254647484e-001;
		xi2[29] = 9.9234295483257617e-001;
		xi2[30] = 9.9854374090973852e-001;
	}



	for (int i = 0; i < points; i++)
	{
		xi[i] = xi[i] * 2 - 1;
	}



	for (int i = 0; i < points2; i++)
	{
		xi1[i] = xi1[i] * 2 - 1;
	}


	for (int i = 0; i < points3; i++)
	{
		xi2[i] = xi2[i] * 2 - 1;
	}



	double *length_vec1 = (double *)malloc(points_length * 3 * sizeof(double));
	double *width_vec1 = (double *)malloc(points_width * 3 * sizeof(double));
	double *height_vec1 = (double *)malloc(points_height * 3 * sizeof(double));
	/*
	if (points_length == 5)
	{
	wi[0] = 0.5688888888888889;
	wi[1] = 0.4786286704993665;
	wi[2] = 0.4786286704993665;
	wi[3] = 0.2369268850561891;
	wi[4] = 0.2369268850561891;

	xi[0] = 0;
	xi[1] = -0.5384693101056831;
	xi[2] = 0.5384693101056831;
	xi[3] = -0.9061798459386640;
	xi[4] = 0.9061798459386640;
	}
	else
	{
	wi[0] = 1;
	wi[1] = 1;

	xi[0] = -0.5773502691896257;
	xi[1] = 0.5773502691896257;
	}

	*/

	for (int i = 0; i < points_length; i++)
	{
		length_vec1[i] = (((fil1->x[1] - fil1->x[0]) / 2)*xi[i])*(fil1->lenvect[0] / (sqrt(pow(fil1->lenvect[0], 2) + pow(fil1->lenvect[1], 2) + pow(fil1->lenvect[2], 2))));
		length_vec1[i + points_length] = (((fil1->y[1] - fil1->y[0]) / 2)*xi[i])*(fil1->lenvect[1] / (sqrt(pow(fil1->lenvect[0], 2) + pow(fil1->lenvect[1], 2) + pow(fil1->lenvect[2], 2))));
		length_vec1[i + points_length * 2] = (((fil1->z[1] - fil1->z[0]) / 2)*xi[i])*(fil1->lenvect[2] / (sqrt(pow(fil1->lenvect[0], 2) + pow(fil1->lenvect[1], 2) + pow(fil1->lenvect[2], 2))));
	}



	/*
	if (points_width == 5)
	{
	wi1[0] = 0.5688888888888889;
	wi1[1] = 0.4786286704993665;
	wi1[2] = 0.4786286704993665;
	wi1[3] = 0.2369268850561891;
	wi1[4] = 0.2369268850561891;

	xi1[0] = 0;
	xi1[1] = -0.5384693101056831;
	xi1[2] = 0.5384693101056831;
	xi1[3] = -0.9061798459386640;
	xi1[4] = 0.9061798459386640;
	}
	else
	{
	wi1[0] = 1;
	wi1[1] = 1;

	xi1[0] = -0.5773502691896257;
	xi1[1] = 0.5773502691896257;
	}
	*/

	for (int i = 0; i < points_width; i++)
	{
		width_vec1[i] = (((fil1->width) / 2)*xi1[i])*width1[0];
		width_vec1[i + points_width] = (((fil1->width) / 2)*xi1[i])*width1[1];
		width_vec1[i + points_width * 2] = (((fil1->width) / 2)*xi1[i])*width1[2];
	}

	/*
	if (points_height == 5)
	{
	wi2[0] = 0.5688888888888889;
	wi2[1] = 0.4786286704993665;
	wi2[2] = 0.4786286704993665;
	wi2[3] = 0.2369268850561891;
	wi2[4] = 0.2369268850561891;

	xi2[0] = 0;
	xi2[1] = -0.5384693101056831;
	xi2[2] = 0.5384693101056831;
	xi2[3] = -0.9061798459386640;
	xi2[4] = 0.9061798459386640;
	}
	else
	{
	wi2[0] = 1;
	wi2[1] = 1;

	xi2[0] = -0.5773502691896257;
	xi2[1] = 0.5773502691896257;
	}
	*/
	for (int i = 0; i < points_height; i++)
	{
		height_vec1[i] = (((fil1->height) / 2)*xi2[i])*height1[0];
		height_vec1[i + points_height] = (((fil1->height) / 2)*xi2[i])*height1[1];
		height_vec1[i + points_height * 2] = (((fil1->height) / 2)*xi2[i])*height1[2];
	}


	//wi2[0] = 0.8888888888888888;
	//wi2[1] = 0.5555555555555556;
	//wi2[2] = 0.5555555555555556;

	//wi2[0] = 1;
	//wi2[1] = 1;





	//xi2[0] = 0;
	//xi2[1] = -0.7745966692414834;
	//xi2[2] = 0.7745966692414834;
	//xi2[0] = -0.5773502691896257;
	//xi2[1] = 0.5773502691896257;








	/*
	width_vec1[i] = (((fil1->width) / 2)*xi[i])*width1[0];
	width_vec1[i + 5] = (((fil1->width) / 2)*xi[i])*width1[1];
	width_vec1[i + 10] = (((fil1->width) / 2)*xi[i])*width1[2];

	height_vec1[i] = (((fil1->height) / 2)*xi[i])*height1[0];
	height_vec1[i + 5] = (((fil1->height) / 2)*xi[i])*height1[1];
	height_vec1[i + 10] = (((fil1->height) / 2)*xi[i])*height1[2];
	*/







	double iii = 0;
	double jjj, kkk;// , lll;
	double x_obs, y_obs, z_obs;
	double dV_R;
	double dV_R_temp;
	//double P0_ij, lplus_ij, R0_ij, Rplus_ij, lmin_ij, Rmin_ij, d_j;
	double d_j;
	//double P0ij_dot_Uij;
	double x_center, y_center, z_center;
	//double x_p, y_p, z_p;
	double l[3];
	double n[3];
	//double p_plus[3];
	//double p_min[3];
	//double p0[3];
	double x_y_z[3];
	double xyz_obs[3];
	double length2[3];
	length2[0] = (fil2->x[1] - fil2->x[0]) / fil2->length;
	length2[1] = (fil2->y[1] - fil2->y[0]) / fil2->length;
	length2[2] = (fil2->z[1] - fil2->z[0]) / fil2->length;
	for (int i = 0; i < points_length; i++)
	{
		jjj = 0;
		for (int j = 0; j < points_width; j++)
		{
			kkk = 0;
			for (int k = 0; k < points_height; k++)
			{
				//lll = 0;
				x_obs = height_vec1[k] + width_vec1[j] + length_vec1[i] + (fil1->x[0] + fil1->x[1]) / 2;
				y_obs = height_vec1[k + points_height] + width_vec1[j + points_width] + length_vec1[i + points_length] + (fil1->y[0] + fil1->y[1]) / 2;
				z_obs = height_vec1[k + points_height * 2] + width_vec1[j + points_width * 2] + length_vec1[i + points_length * 2] + (fil1->z[0] + fil1->z[1]) / 2;
				//Start Integral



				//start first edge

				xyz_obs[0] = x_obs;
				xyz_obs[1] = y_obs;
				xyz_obs[2] = z_obs;
				xyz_point(x_y_z, length2, xyz_obs, fil2->x[0] + height2[0], fil2->y[0] + height2[1], fil2->z[0] + height2[2]);


				x_center = fil2->x[0] - height2[0] * fil2->height / 2;
				y_center = fil2->y[0] - height2[1] * fil2->height / 2;
				z_center = fil2->z[0] - height2[2] * fil2->height / 2;


				l[0] = -width2[0];
				l[1] = -width2[1];
				l[2] = -width2[2];


				n[0] = (fil2->x[0] - fil2->x[1]) / fil2->length;
				n[1] = (fil2->y[0] - fil2->y[1]) / fil2->length;
				n[2] = (fil2->z[0] - fil2->z[1]) / fil2->length;



				dV_R_temp = dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, -fil2->width*width2[0] / 2, -fil2->width*width2[1] / 2, -fil2->width*width2[2] / 2);

				//end first edge

				//start second edge


				x_center = fil2->x[0] + height2[0] * fil2->height / 2;//change from first +/-
				y_center = fil2->y[0] + height2[1] * fil2->height / 2;
				z_center = fil2->z[0] + height2[2] * fil2->height / 2;



				//double P0 = sqrt(pow(x_proj - x_p,2)+pow(y_proj - y_p,2)+pow(z_proj - z_p,2));



				l[0] = width2[0];
				l[1] = width2[1];
				l[2] = width2[2];




				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, fil2->width*width2[0] / 2, fil2->width*width2[1] / 2, fil2->width*width2[2] / 2);

				//end second edge

				//start third edge

				x_center = fil2->x[0] - width2[0] * fil2->width / 2;//change height->width
				y_center = fil2->y[0] - width2[1] * fil2->width / 2;
				z_center = fil2->z[0] - width2[2] * fil2->width / 2;




				//double P0 = sqrt(pow(x_proj - x_p,2)+pow(y_proj - y_p,2)+pow(z_proj - z_p,2));



				l[0] = height2[0];
				l[1] = height2[1];
				l[2] = height2[2];


				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, fil2->height*height2[0] / 2, fil2->height*height2[1] / 2, fil2->height*height2[2] / 2);

				//end third edge



				//start forth edge


				x_center = fil2->x[0] + width2[0] * fil2->width / 2;//change height->width
				y_center = fil2->y[0] + width2[1] * fil2->width / 2;
				z_center = fil2->z[0] + width2[2] * fil2->width / 2;


				//double P0 = sqrt(pow(x_proj - x_p,2)+pow(y_proj - y_p,2)+pow(z_proj - z_p,2));


				l[0] = -height2[0];
				l[1] = -height2[1];
				l[2] = -height2[2];





				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, -fil2->height*height2[0] / 2, -fil2->height*height2[1] / 2, -fil2->height*height2[2] / 2);



				//end forth edge					
				d_j = (x_obs - x_center)*n[0] + (y_obs - y_center)*n[1] + (z_obs - z_center)*n[2];
				dV_R_temp = d_j*dV_R_temp;
				dV_R = dV_R_temp;

				//start fifth edge
				n[0] = (fil2->x[1] - fil2->x[0]) / fil2->length;
				n[1] = (fil2->y[1] - fil2->y[0]) / fil2->length;
				n[2] = (fil2->z[1] - fil2->z[0]) / fil2->length;



				//xyz_obs[0] = x_obs;
				//xyz_obs[1] = y_obs;
				//xyz_obs[2] = z_obs;
				xyz_point(x_y_z, length2, xyz_obs, fil2->x[1], fil2->y[1], fil2->z[1]);

				x_center = fil2->x[1] - height2[0] * fil2->height / 2;
				y_center = fil2->y[1] - height2[1] * fil2->height / 2;
				z_center = fil2->z[1] - height2[2] * fil2->height / 2;


				l[0] = width2[0];
				l[1] = width2[1];
				l[2] = width2[2];

				dV_R_temp = dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, fil2->width*width2[0] / 2, fil2->width*width2[1] / 2, fil2->width*width2[2] / 2);
				//end fifth edge

				//start sixth edge


				x_center = fil2->x[1] + height2[0] * fil2->height / 2;//change from first +/-
				y_center = fil2->y[1] + height2[1] * fil2->height / 2;
				z_center = fil2->z[1] + height2[2] * fil2->height / 2;

				//double P0 = sqrt(pow(x_proj - x_p,2)+pow(y_proj - y_p,2)+pow(z_proj - z_p,2));


				l[0] = -width2[0];
				l[1] = -width2[1];
				l[2] = -width2[2];


				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, -fil2->width*width2[0] / 2, -fil2->width*width2[1] / 2, -fil2->width*width2[2] / 2);


				//end sixth edge

				//start seventh edge


				x_center = fil2->x[1] - width2[0] * fil2->width / 2;//change height->width
				y_center = fil2->y[1] - width2[1] * fil2->width / 2;
				z_center = fil2->z[1] - width2[2] * fil2->width / 2;



				//double P0 = sqrt(pow(x_proj - x_p,2)+pow(y_proj - y_p,2)+pow(z_proj - z_p,2));


				l[0] = -height2[0];
				l[1] = -height2[1];
				l[2] = -height2[2];


				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, -fil2->height*height2[0] / 2, -fil2->height*height2[1] / 2, -fil2->height*height2[2] / 2);

				//end seventh edge

				//start eighth edge

				x_center = fil2->x[1] + width2[0] * fil2->width / 2;//change height->width
				y_center = fil2->y[1] + width2[1] * fil2->width / 2;
				z_center = fil2->z[1] + width2[2] * fil2->width / 2;




				//double P0 = sqrt(pow(x_proj - x_p,2)+pow(y_proj - y_p,2)+pow(z_proj - z_p,2));



				l[0] = height2[0];
				l[1] = height2[1];
				l[2] = height2[2];



				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, fil2->height*height2[0] / 2, fil2->height*height2[1] / 2, fil2->height*height2[2] / 2);
				//end eigth edge
				d_j = (x_obs - x_center)*n[0] + (y_obs - y_center)*n[1] + (z_obs - z_center)*n[2];
				dV_R_temp = d_j*dV_R_temp;
				dV_R += dV_R_temp;



				//start ninth edge

				x_center = (fil2->x[0] + fil2->x[1]) / 2 + width2[0] * fil2->width / 2 - height2[0] * fil2->height / 2;
				y_center = (fil2->y[0] + fil2->y[1]) / 2 + width2[1] * fil2->width / 2 - height2[1] * fil2->height / 2;
				z_center = (fil2->z[0] + fil2->z[1]) / 2 + width2[2] * fil2->width / 2 - height2[2] * fil2->height / 2;

				xyz_obs[0] = x_obs;
				xyz_obs[1] = y_obs;
				xyz_obs[2] = z_obs;
				xyz_point(x_y_z, width2, xyz_obs, x_center, y_center, z_center);

				n[0] = width2[0];
				n[1] = width2[1];
				n[2] = width2[2];


				l[0] = (fil2->x[0] - fil2->x[1]) / fil2->length;
				l[1] = (fil2->y[0] - fil2->y[1]) / fil2->length;
				l[2] = (fil2->z[0] - fil2->z[1]) / fil2->length;

				dV_R_temp = dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, -fil2->length*length2[0] / 2, -fil2->length*length2[1] / 2, -fil2->length*length2[2] / 2);

				//end ninth edge

				//start tenth edge

				x_center = (fil2->x[0] + fil2->x[1]) / 2 + width2[0] * fil2->width / 2 + height2[0] * fil2->height / 2;
				y_center = (fil2->y[0] + fil2->y[1]) / 2 + width2[1] * fil2->width / 2 + height2[1] * fil2->height / 2;
				z_center = (fil2->z[0] + fil2->z[1]) / 2 + width2[2] * fil2->width / 2 + height2[2] * fil2->height / 2;




				//l[0] = fil2->lenvect[0];
				//l[1] = fil2->lenvect[1];
				//l[2] = fil2->lenvect[2];

				l[0] = (fil2->x[1] - fil2->x[0]) / fil2->length;
				l[1] = (fil2->y[1] - fil2->y[0]) / fil2->length;
				l[2] = (fil2->z[1] - fil2->z[0]) / fil2->length;

				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, fil2->length*length2[0] / 2, fil2->length*length2[1] / 2, fil2->length*length2[2] / 2);
				//end tenth edge


				//start eleventh edge


				x_center = fil2->x[0] + width2[0] * fil2->width / 2;
				y_center = fil2->y[0] + width2[1] * fil2->width / 2;
				z_center = fil2->z[0] + width2[2] * fil2->width / 2;




				l[0] = height2[0];
				l[1] = height2[1];
				l[2] = height2[2];


				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, fil2->height*height2[0] / 2, fil2->height*height2[1] / 2, fil2->height*height2[2] / 2);

				//end eleventh edge


				//start twelfth edge


				x_center = fil2->x[1] + width2[0] * fil2->width / 2;
				y_center = fil2->y[1] + width2[1] * fil2->width / 2;
				z_center = fil2->z[1] + width2[2] * fil2->width / 2;




				l[0] = -height2[0];
				l[1] = -height2[1];
				l[2] = -height2[2];


				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, -fil2->height*height2[0] / 2, -fil2->height*height2[1] / 2, -fil2->height*height2[2] / 2);
				//end twelth edge
				d_j = (x_obs - x_center)*n[0] + (y_obs - y_center)*n[1] + (z_obs - z_center)*n[2];
				dV_R_temp = d_j*dV_R_temp;
				dV_R += dV_R_temp;




				//start thirtenth edge

				x_center = (fil2->x[0] + fil2->x[1]) / 2 - width2[0] * fil2->width / 2 - height2[0] * fil2->height / 2;
				y_center = (fil2->y[0] + fil2->y[1]) / 2 - width2[1] * fil2->width / 2 - height2[1] * fil2->height / 2;
				z_center = (fil2->z[0] + fil2->z[1]) / 2 - width2[2] * fil2->width / 2 - height2[2] * fil2->height / 2;

				//xyz_obs[0] = x_obs;
				//xyz_obs[1] = y_obs;
				//xyz_obs[2] = z_obs;
				xyz_point(x_y_z, width2, xyz_obs, x_center, y_center, z_center);

				n[0] = -width2[0];
				n[1] = -width2[1];
				n[2] = -width2[2];

				//l[0] = fil2->lenvect[0];
				//l[1] = fil2->lenvect[1];
				//l[2] = fil2->lenvect[2];

				l[0] = (fil2->x[1] - fil2->x[0]) / fil2->length;
				l[1] = (fil2->y[1] - fil2->y[0]) / fil2->length;
				l[2] = (fil2->z[1] - fil2->z[0]) / fil2->length;

				dV_R_temp = dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, fil2->length*length2[0] / 2, fil2->length*length2[1] / 2, fil2->length*length2[2] / 2);

				//end thirtenth edge

				//start fourteenth edge


				x_center = (fil2->x[0] + fil2->x[1]) / 2 - width2[0] * fil2->width / 2 + height2[0] * fil2->height / 2;
				y_center = (fil2->y[0] + fil2->y[1]) / 2 - width2[1] * fil2->width / 2 + height2[1] * fil2->height / 2;
				z_center = (fil2->z[0] + fil2->z[1]) / 2 - width2[2] * fil2->width / 2 + height2[2] * fil2->height / 2;



				//l[0] = -fil2->lenvect[0];
				//l[1] = -fil2->lenvect[1];
				//l[2] = -fil2->lenvect[2];

				l[0] = (fil2->x[0] - fil2->x[1]) / fil2->length;
				l[1] = (fil2->y[0] - fil2->y[1]) / fil2->length;
				l[2] = (fil2->z[0] - fil2->z[1]) / fil2->length;

				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, -fil2->length*length2[0] / 2, -fil2->length*length2[1] / 2, -fil2->length*length2[2] / 2);
				//end fourteenth edge

				//start fifteenth edge

				x_center = fil2->x[0] - width2[0] * fil2->width / 2;
				y_center = fil2->y[0] - width2[1] * fil2->width / 2;
				z_center = fil2->z[0] - width2[2] * fil2->width / 2;




				l[0] = -height2[0];
				l[1] = -height2[1];
				l[2] = -height2[2];

				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, -fil2->height*height2[0] / 2, -fil2->height*height2[1] / 2, -fil2->height*height2[2] / 2);
				//end fifteenth edge

				//start sixteenth edge

				x_center = fil2->x[1] - width2[0] * fil2->width / 2;
				y_center = fil2->y[1] - width2[1] * fil2->width / 2;
				z_center = fil2->z[1] - width2[2] * fil2->width / 2;



				l[0] = height2[0];
				l[1] = height2[1];
				l[2] = height2[2];

				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, fil2->height*height2[0] / 2, fil2->height*height2[1] / 2, fil2->height*height2[2] / 2);
				//end sixteenth edge
				d_j = (x_obs - x_center)*n[0] + (y_obs - y_center)*n[1] + (z_obs - z_center)*n[2];
				dV_R_temp = d_j*dV_R_temp;
				dV_R += dV_R_temp;



				//start seventeenth edge



				x_center = fil2->x[0] - height2[0] * fil2->height / 2;
				y_center = fil2->y[0] - height2[1] * fil2->height / 2;
				z_center = fil2->z[0] - height2[2] * fil2->height / 2;

				xyz_obs[0] = x_obs;
				xyz_obs[1] = y_obs;
				xyz_obs[2] = z_obs;
				xyz_point(x_y_z, height2, xyz_obs, x_center, y_center, z_center);



				n[0] = -height2[0];
				n[1] = -height2[1];
				n[2] = -height2[2];

				l[0] = width2[0];
				l[1] = width2[1];
				l[2] = width2[2];

				dV_R_temp = dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, fil2->width*width2[0] / 2, fil2->width*width2[1] / 2, fil2->width*width2[2] / 2);

				//end seventeenth edge

				//start eighteenth edge

				x_center = fil2->x[1] - height2[0] * fil2->height / 2;
				y_center = fil2->y[1] - height2[1] * fil2->height / 2;
				z_center = fil2->z[1] - height2[2] * fil2->height / 2;




				l[0] = -width2[0];
				l[1] = -width2[1];
				l[2] = -width2[2];

				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, -fil2->width*width2[0] / 2, -fil2->width*width2[1] / 2, -fil2->width*width2[2] / 2);
				//end eighteenth edge

				//start nineteenth edge

				x_center = (fil2->x[0] + fil2->x[1]) / 2 - height2[0] * fil2->height / 2 + width2[0] * fil2->width / 2;
				y_center = (fil2->y[0] + fil2->y[1]) / 2 - height2[1] * fil2->height / 2 + width2[1] * fil2->width / 2;
				z_center = (fil2->z[0] + fil2->z[1]) / 2 - height2[2] * fil2->height / 2 + width2[2] * fil2->width / 2;



				//l[0] = fil2->lenvect[0];
				//l[1] = fil2->lenvect[1];
				//l[2] = fil2->lenvect[2];

				l[0] = (fil2->x[1] - fil2->x[0]) / fil2->length;
				l[1] = (fil2->y[1] - fil2->y[0]) / fil2->length;
				l[2] = (fil2->z[1] - fil2->z[0]) / fil2->length;

				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, fil2->length*length2[0] / 2, fil2->length*length2[1] / 2, fil2->length*length2[2] / 2);
				//end nineteenth edge

				//start twentieth edge


				x_center = (fil2->x[0] + fil2->x[1]) / 2 - height2[0] * fil2->height / 2 - width2[0] * fil2->width / 2;
				y_center = (fil2->y[0] + fil2->y[1]) / 2 - height2[1] * fil2->height / 2 - width2[1] * fil2->width / 2;
				z_center = (fil2->z[0] + fil2->z[1]) / 2 - height2[2] * fil2->height / 2 - width2[2] * fil2->width / 2;





				//l[0] = -fil2->lenvect[0];
				//l[1] = -fil2->lenvect[1];
				//l[2] = -fil2->lenvect[2];

				l[0] = (fil2->x[0] - fil2->x[1]) / fil2->length;
				l[1] = (fil2->y[0] - fil2->y[1]) / fil2->length;
				l[2] = (fil2->z[0] - fil2->z[1]) / fil2->length;

				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, -fil2->length*length2[0] / 2, -fil2->length*length2[1] / 2, -fil2->length*length2[2] / 2);
				//end twentieth edge
				d_j = (x_obs - x_center)*n[0] + (y_obs - y_center)*n[1] + (z_obs - z_center)*n[2];
				dV_R_temp = d_j*dV_R_temp;
				dV_R += dV_R_temp;



				//start twentiefirst edge

				x_center = fil2->x[0] + height2[0] * fil2->height / 2;
				y_center = fil2->y[0] + height2[1] * fil2->height / 2;
				z_center = fil2->z[0] + height2[2] * fil2->height / 2;

				xyz_obs[0] = x_obs;
				xyz_obs[1] = y_obs;
				xyz_obs[2] = z_obs;
				xyz_point(x_y_z, height2, xyz_obs, x_center, y_center, z_center);

				n[0] = height2[0];
				n[1] = height2[1];
				n[2] = height2[2];

				l[0] = -width2[0];
				l[1] = -width2[1];
				l[2] = -width2[2];

				dV_R_temp = dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, -fil2->width*width2[0] / 2, -fil2->width*width2[1] / 2, -fil2->width*width2[2] / 2);

				// end twentiefirst edge


				//start twentiesecond edge

				x_center = fil2->x[1] + height2[0] * fil2->height / 2;
				y_center = fil2->y[1] + height2[1] * fil2->height / 2;
				z_center = fil2->z[1] + height2[2] * fil2->height / 2;



				l[0] = width2[0];
				l[1] = width2[1];
				l[2] = width2[2];

				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, fil2->width*width2[0] / 2, fil2->width*width2[1] / 2, fil2->width*width2[2] / 2);
				//end twentiesecond edge

				//start twentiethird edge

				x_center = (fil2->x[0] + fil2->x[1]) / 2 + height2[0] * fil2->height / 2 - width2[0] * fil2->width / 2;
				y_center = (fil2->y[0] + fil2->y[1]) / 2 + height2[1] * fil2->height / 2 - width2[1] * fil2->width / 2;
				z_center = (fil2->z[0] + fil2->z[1]) / 2 + height2[2] * fil2->height / 2 - width2[2] * fil2->width / 2;



				//l[0] = fil2->lenvect[0];
				//l[1] = fil2->lenvect[1];
				//l[2] = fil2->lenvect[2];

				l[0] = (fil2->x[1] - fil2->x[0]) / fil2->length;
				l[1] = (fil2->y[1] - fil2->y[0]) / fil2->length;
				l[2] = (fil2->z[1] - fil2->z[0]) / fil2->length;

				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, fil2->length*length2[0] / 2, fil2->length*length2[1] / 2, fil2->length*length2[2] / 2);
				//end twentiethird edge


				//start twentieforth edge

				x_center = (fil2->x[0] + fil2->x[1]) / 2 + height2[0] * fil2->height / 2 + width2[0] * fil2->width / 2;
				y_center = (fil2->y[0] + fil2->y[1]) / 2 + height2[1] * fil2->height / 2 + width2[1] * fil2->width / 2;
				z_center = (fil2->z[0] + fil2->z[1]) / 2 + height2[2] * fil2->height / 2 + width2[2] * fil2->width / 2;




				//l[0] = -fil2->lenvect[0];
				//l[1] = -fil2->lenvect[1];
				//l[2] = -fil2->lenvect[2];

				l[0] = (fil2->x[0] - fil2->x[1]) / fil2->length;
				l[1] = (fil2->y[0] - fil2->y[1]) / fil2->length;
				l[2] = (fil2->z[0] - fil2->z[1]) / fil2->length;

				dV_R_temp += dV_R_temp1(x_y_z[0], x_y_z[1], x_y_z[2], x_center, y_center, z_center, x_obs, y_obs, z_obs, l, n
					, -fil2->length*length2[0] / 2, -fil2->length*length2[1] / 2, -fil2->length*length2[2] / 2);
				//end twentieforth edge
				d_j = (x_obs - x_center)*n[0] + (y_obs - y_center)*n[1] + (z_obs - z_center)*n[2];
				dV_R_temp = d_j*dV_R_temp;
				dV_R += dV_R_temp;
				//End Integral
				//lll = 0.5*dV_R;
				kkk += wi2[k] * 0.5*dV_R;
			}
			jjj += wi1[j] * kkk;
		}
		iii += wi[i] * jjj;
	}


	//double tttttt = (MU0 / (4 * PI *fil1->area*fil2->area)) * a_b * iii;
	//if (tttttt == 0)
	//{
	//	int dsjhfdsjf = 0;
	//}
	//return (MU0 / (4 * PI */**fil1->area*/fil2->area)) * a_b * iii;

	free(wi);
	free(wi1);
	free(wi2);
	free(xi);
	free(xi1);
	free(xi2);
	free(length_vec1);
	free(width_vec1);
	free(height_vec1);

	return   iii / (fil2->area*fil2->length * 8);
}


//From here is adapted from mutual.c until  mutualfill3
int realcos_error = 0;

double magdiff2(FILAMENT *fil1, int node1, FILAMENT *fil2, int node2) {
	return (SQUARE(fil1->x[node1] - fil2->x[node2])
		+ SQUARE(fil1->y[node1] - fil2->y[node2])
		+ SQUARE(fil1->z[node1] - fil2->z[node2]));
}

/* this gives the mutual inductance of two filaments who represent */
/* opposite sides of a rectangle */

/*
double mut_rect(double len, double d) {
	double temp, temp1;

	temp = sqrt(len * len + d * d);
	temp1 = len * asinh(len / d);
	return temp - temp1;
}
*/



/*Similar to Aanalitic but now Gaussian integration is used on the inner and outer 
integral this is typically used for well separated filaments.*/
double mutual_quadrature(FILAMENT *fil1, FILAMENT *fil2, int points, int points2, int points3, int points4, int points5, int points6)//this asumes orientation in only one axis 
{
	

	int points_length = points;
	int points_width = points2;
	int points_height = points3;


	double width1[3];
	double width2[3];
	double height1[3];
	double height2[3];
	get_wid(fil1, width1);
	get_wid(fil2, width2);

	get_height(fil1, width1, height1);
	get_height(fil2, width2, height2);

	double *wi = (double *)malloc(points_length * sizeof(double));
	double *wi1 = (double *)malloc(points_width * sizeof(double));
	double *wi2 = (double *)malloc(points_height * sizeof(double));

	double *xi = (double *)malloc(points_length * sizeof(double));
	double *xi1 = (double *)malloc(points_width * sizeof(double));
	double *xi2 = (double *)malloc(points_height * sizeof(double));

	double *length_vec1 = (double *)malloc(points_length * 3 * sizeof(double));
	double *width_vec1 = (double *)malloc(points_width * 3 * sizeof(double));
	double *height_vec1 = (double *)malloc(points_height * 3 * sizeof(double));



	if (points == 1)
	{
		wi[0] = 1.0000000000000000e+000;
	}
	else if (points == 2)
	{
		wi[0] = 4.9999999999999989e-001;
		wi[1] = 4.9999999999999989e-001;
	}
	else if (points == 3)
	{
		wi[0] = 2.7777777777777768e-001;
		wi[1] = 4.4444444444444425e-001;
		wi[2] = 2.7777777777777779e-001;

	}
	else if (points == 4)
	{
		wi[0] = 1.7392742256872706e-001;
		wi[1] = 3.2607257743127316e-001;
		wi[2] = 3.2607257743127299e-001;
		wi[3] = 1.7392742256872701e-001;
	}
	else if (points == 5)
	{
		wi[0] = 1.1846344252809456e-001;
		wi[1] = 2.3931433524968337e-001;
		wi[2] = 2.8444444444444444e-001;
		wi[3] = 2.3931433524968343e-001;
		wi[4] = 1.1846344252809445e-001;
	}
	else if (points == 6)
	{
		wi[0] = 8.5662246189585289e-002;
		wi[1] = 1.8038078652406928e-001;
		wi[2] = 2.3395696728634516e-001;
		wi[3] = 2.3395696728634591e-001;
		wi[4] = 1.8038078652406928e-001;
		wi[5] = 8.5662246189585220e-002;
	}
	else if (points == 7)
	{
		wi[0] = 6.4742483084434824e-002;
		wi[1] = 1.3985269574463838e-001;
		wi[2] = 1.9091502525255943e-001;
		wi[3] = 2.0897959183673431e-001;
		wi[4] = 1.9091502525255952e-001;
		wi[5] = 1.3985269574463827e-001;
		wi[6] = 6.4742483084434907e-002;
	}



	if (points2 == 1)
	{
		wi1[0] = 1.0000000000000000e+000;
	}
	else if (points2 == 2)
	{
		wi1[0] = 4.9999999999999989e-001;
		wi1[1] = 4.9999999999999989e-001;
	}
	else if (points2 == 3)
	{
		wi1[0] = 2.7777777777777768e-001;
		wi1[1] = 4.4444444444444425e-001;
		wi1[2] = 2.7777777777777779e-001;

	}
	else if (points2 == 4)
	{
		wi1[0] = 1.7392742256872706e-001;
		wi1[1] = 3.2607257743127316e-001;
		wi1[2] = 3.2607257743127299e-001;
		wi1[3] = 1.7392742256872701e-001;
	}
	else if (points2 == 5)
	{
		wi1[0] = 1.1846344252809456e-001;
		wi1[1] = 2.3931433524968337e-001;
		wi1[2] = 2.8444444444444444e-001;
		wi1[3] = 2.3931433524968343e-001;
		wi1[4] = 1.1846344252809445e-001;
	}
	else if (points2 == 6)
	{
		wi1[0] = 8.5662246189585289e-002;
		wi1[1] = 1.8038078652406928e-001;
		wi1[2] = 2.3395696728634516e-001;
		wi1[3] = 2.3395696728634591e-001;
		wi1[4] = 1.8038078652406928e-001;
		wi1[5] = 8.5662246189585220e-002;
	}
	else if (points2 == 7)
	{

		wi1[0] = 6.4742483084434824e-002;
		wi1[1] = 1.3985269574463838e-001;
		wi1[2] = 1.9091502525255943e-001;
		wi1[3] = 2.0897959183673431e-001;
		wi1[4] = 1.9091502525255952e-001;
		wi1[5] = 1.3985269574463827e-001;
		wi1[6] = 6.4742483084434907e-002;

	}




	if (points3 == 1)
	{
		wi2[0] = 1.0000000000000000e+000;
	}
	else if (points3 == 2)
	{
		wi2[0] = 4.9999999999999989e-001;
		wi2[1] = 4.9999999999999989e-001;
	}
	else if (points3 == 3)
	{
		wi2[0] = 2.7777777777777768e-001;
		wi2[1] = 4.4444444444444425e-001;
		wi2[2] = 2.7777777777777779e-001;
	}
	else if (points3 == 4)
	{
		wi2[0] = 1.7392742256872706e-001;
		wi2[1] = 3.2607257743127316e-001;
		wi2[2] = 3.2607257743127299e-001;
		wi2[3] = 1.7392742256872701e-001;
	}
	else if (points3 == 5)
	{
		wi2[0] = 1.1846344252809456e-001;
		wi2[1] = 2.3931433524968337e-001;
		wi2[2] = 2.8444444444444444e-001;
		wi2[3] = 2.3931433524968343e-001;
		wi2[4] = 1.1846344252809445e-001;
	}
	else if (points3 == 6)
	{
		wi2[0] = 8.5662246189585289e-002;
		wi2[1] = 1.8038078652406928e-001;
		wi2[2] = 2.3395696728634516e-001;
		wi2[3] = 2.3395696728634591e-001;
		wi2[4] = 1.8038078652406928e-001;
		wi2[5] = 8.5662246189585220e-002;
	}
	else if (points3 == 7)
	{

		wi2[0] = 6.4742483084434824e-002;
		wi2[1] = 1.3985269574463838e-001;
		wi2[2] = 1.9091502525255943e-001;
		wi2[3] = 2.0897959183673431e-001;
		wi2[4] = 1.9091502525255952e-001;
		wi2[5] = 1.3985269574463827e-001;
		wi2[6] = 6.4742483084434907e-002;
	}



	if (points == 1)
	{
		xi[0] = 5.0000000000000000e-001;
	}
	else if (points == 2)
	{
		xi[0] = 2.1132486540518708e-001;
		xi[1] = 7.8867513459481287e-001;
	}
	else if (points == 3)
	{
		xi[0] = 1.1270166537925824e-001;
		xi[1] = 5.0000000000000000e-001;
		xi[2] = 8.8729833462074170e-001;

	}
	else if (points == 4)
	{
		xi[0] = 6.9431844202973714e-002;
		xi[1] = 3.3000947820757182e-001;
		xi[2] = 6.6999052179242813e-001;
		xi[3] = 9.3056815579702634e-001;
	}
	else if (points == 5)
	{
		xi[0] = 4.6910077030668074e-002;
		xi[1] = 2.3076534494715845e-001;
		xi[2] = 5.0000000000000000e-001;
		xi[3] = 7.6923465505284150e-001;
		xi[4] = 9.5308992296933193e-001;
	}
	else if (points == 6)
	{
		xi[0] = 3.3765242898423975e-002;
		xi[1] = 1.6939530676686776e-001;
		xi[2] = 3.8069040695840151e-001;
		xi[3] = 6.1930959304159849e-001;
		xi[4] = 8.3060469323313224e-001;
		xi[5] = 9.6623475710157614e-001;
	}
	else if (points == 7)
	{

		xi[0] = 2.5446043828620812e-002;
		xi[1] = 1.2923440720030277e-001;
		xi[2] = 2.9707742431130146e-001;
		xi[3] = 4.9999999999999994e-001;
		xi[4] = 7.0292257568869854e-001;
		xi[5] = 8.7076559279969734e-001;
		xi[6] = 9.7455395617137919e-001;
	}




	if (points2 == 1)
	{
		xi1[0] = 5.0000000000000000e-001;
	}
	else if (points2 == 2)
	{
		xi1[0] = 2.1132486540518708e-001;
		xi1[1] = 7.8867513459481287e-001;
	}
	else if (points2 == 3)
	{
		xi1[0] = 1.1270166537925824e-001;
		xi1[1] = 5.0000000000000000e-001;
		xi1[2] = 8.8729833462074170e-001;

	}
	else if (points2 == 4)
	{
		xi1[0] = 6.9431844202973714e-002;
		xi1[1] = 3.3000947820757182e-001;
		xi1[2] = 6.6999052179242813e-001;
		xi1[3] = 9.3056815579702634e-001;
	}
	else if (points2 == 5)
	{
		xi1[0] = 4.6910077030668074e-002;
		xi1[1] = 2.3076534494715845e-001;
		xi1[2] = 5.0000000000000000e-001;
		xi1[3] = 7.6923465505284150e-001;
		xi1[4] = 9.5308992296933193e-001;
	}
	else if (points2 == 6)
	{
		xi1[0] = 3.3765242898423975e-002;
		xi1[1] = 1.6939530676686776e-001;
		xi1[2] = 3.8069040695840151e-001;
		xi1[3] = 6.1930959304159849e-001;
		xi1[4] = 8.3060469323313224e-001;
		xi1[5] = 9.6623475710157614e-001;
	}
	else if (points2 == 7)
	{
		xi1[0] = 2.5446043828620812e-002;
		xi1[1] = 1.2923440720030277e-001;
		xi1[2] = 2.9707742431130146e-001;
		xi1[3] = 4.9999999999999994e-001;
		xi1[4] = 7.0292257568869854e-001;
		xi1[5] = 8.7076559279969734e-001;
		xi1[6] = 9.7455395617137919e-001;
	}




	if (points3 == 1)
	{
		xi2[0] = 5.0000000000000000e-001;
	}
	else if (points3 == 2)
	{
		xi2[0] = 2.1132486540518708e-001;
		xi2[1] = 7.8867513459481287e-001;
	}
	else if (points3 == 3)
	{
		xi2[0] = 1.1270166537925824e-001;
		xi2[1] = 5.0000000000000000e-001;
		xi2[2] = 8.8729833462074170e-001;
	}
	else if (points3 == 4)
	{
		xi2[0] = 6.9431844202973714e-002;
		xi2[1] = 3.3000947820757182e-001;
		xi2[2] = 6.6999052179242813e-001;
		xi2[3] = 9.3056815579702634e-001;
	}
	else if (points3 == 5)
	{
		xi2[0] = 4.6910077030668074e-002;
		xi2[1] = 2.3076534494715845e-001;
		xi2[2] = 5.0000000000000000e-001;
		xi2[3] = 7.6923465505284150e-001;
		xi2[4] = 9.5308992296933193e-001;
	}
	else if (points3 == 6)
	{
		xi2[0] = 3.3765242898423975e-002;
		xi2[1] = 1.6939530676686776e-001;
		xi2[2] = 3.8069040695840151e-001;
		xi2[3] = 6.1930959304159849e-001;
		xi2[4] = 8.3060469323313224e-001;
		xi2[5] = 9.6623475710157614e-001;
	}
	else if (points3 == 7)
	{
		xi2[0] = 2.5446043828620812e-002;
		xi2[1] = 1.2923440720030277e-001;
		xi2[2] = 2.9707742431130146e-001;
		xi2[3] = 4.9999999999999994e-001;
		xi2[4] = 7.0292257568869854e-001;
		xi2[5] = 8.7076559279969734e-001;
		xi2[6] = 9.7455395617137919e-001;
	}


	//Next 6 for loops are redundant. It is recommended to change Gaussian quadrature values.
	for (int i = 0; i < points; i++)
	{
		wi[i] *= 2;
	}


	for (int i = 0; i < points2; i++)
	{
		wi1[i] *= 2;
	}


	for (int i = 0; i < points3; i++)
	{
		wi2[i] *= 2;

	}

	for (int i = 0; i < points; i++)
	{
		xi[i] = xi[i] * 2 - 1;
	}

	for (int i = 0; i < points2; i++)
	{
		xi1[i] = xi1[i] * 2 - 1;
	}


	for (int i = 0; i < points3; i++)
	{
		xi2[i] = xi2[i] * 2 - 1;
	}

	/*
	if (points_length == 5)
	{
	wi[0] = 0.5688888888888889;
	wi[1] = 0.4786286704993665;
	wi[2] = 0.4786286704993665;
	wi[3] = 0.2369268850561891;
	wi[4] = 0.2369268850561891;

	xi[0] = 0;
	xi[1] = -0.5384693101056831;
	xi[2] = 0.5384693101056831;
	xi[3] = -0.9061798459386640;
	xi[4] = 0.9061798459386640;
	}
	else
	{
	wi[0] = 1;
	wi[1] = 1;

	xi[0] = -0.5773502691896257;
	xi[1] = 0.5773502691896257;
	}

	*/




	for (int i = 0; i < points_length; i++)
	{
		length_vec1[i] = (((fil1->x[1] - fil1->x[0]) / 2)*xi[i])*(fil1->lenvect[0] / (sqrt(pow(fil1->lenvect[0], 2) + pow(fil1->lenvect[1], 2) + pow(fil1->lenvect[2], 2))));
		length_vec1[i + points_length] = (((fil1->y[1] - fil1->y[0]) / 2)*xi[i])*(fil1->lenvect[1] / (sqrt(pow(fil1->lenvect[0], 2) + pow(fil1->lenvect[1], 2) + pow(fil1->lenvect[2], 2))));
		length_vec1[i + points_length * 2] = (((fil1->z[1] - fil1->z[0]) / 2)*xi[i])*(fil1->lenvect[2] / (sqrt(pow(fil1->lenvect[0], 2) + pow(fil1->lenvect[1], 2) + pow(fil1->lenvect[2], 2))));
	}


	/*

	if (points_width == 5)
	{
	wi1[0] = 0.5688888888888889;
	wi1[1] = 0.4786286704993665;
	wi1[2] = 0.4786286704993665;
	wi1[3] = 0.2369268850561891;
	wi1[4] = 0.2369268850561891;

	xi1[0] = 0;
	xi1[1] = -0.5384693101056831;
	xi1[2] = 0.5384693101056831;
	xi1[3] = -0.9061798459386640;
	xi1[4] = 0.9061798459386640;
	}
	else
	{
	wi1[0] = 1;
	wi1[1] = 1;

	xi1[0] = -0.5773502691896257;
	xi1[1] = 0.5773502691896257;
	}
	*/


	for (int i = 0; i < points_width; i++)
	{
		width_vec1[i] = (((fil1->width) / 2)*xi1[i])*width1[0];
		width_vec1[i + points_width] = (((fil1->width) / 2)*xi1[i])*width1[1];
		width_vec1[i + points_width * 2] = (((fil1->width) / 2)*xi1[i])*width1[2];
	}

	/*
	if (points_height == 5)
	{
	wi2[0] = 0.5688888888888889;
	wi2[1] = 0.4786286704993665;
	wi2[2] = 0.4786286704993665;
	wi2[3] = 0.2369268850561891;
	wi2[4] = 0.2369268850561891;

	xi2[0] = 0;
	xi2[1] = -0.5384693101056831;
	xi2[2] = 0.5384693101056831;
	xi2[3] = -0.9061798459386640;
	xi2[4] = 0.9061798459386640;
	}
	else
	{
	wi2[0] = 1;
	wi2[1] = 1;

	xi2[0] = -0.5773502691896257;
	xi2[1] = 0.5773502691896257;
	}
	*/
	for (int i = 0; i < points_height; i++)
	{
		height_vec1[i] = (((fil1->height) / 2)*xi2[i])*height1[0];
		height_vec1[i + points_height] = (((fil1->height) / 2)*xi2[i])*height1[1];
		height_vec1[i + points_height * 2] = (((fil1->height) / 2)*xi2[i])*height1[2];

	}







	//double length_vec2[15];
	//double width_vec2[15];
	//double height_vec2[15];





	int points_length2 = points4;
	int points_width2 = points5;
	int points_height2 = points6;

	/*
	if (fil2->length > fil2->width && fil2->length > fil2->height)
	{
	points_length2 = 5;
	if (fil2->width * 2 > fil2->length )
	{
	points_width2 = 5;
	}
	else
	{
	points_width2 = 2;
	}
	if (fil2->height * 2 > fil2->length)
	{
	points_height2 = 5;
	}
	else
	{
	points_height2 = 2;
	}

	}

	else if (fil2->width > fil2->height)
	{
	points_width2 = 5;

	if (fil2->length * 2 > fil2->width)
	{
	points_length2 = 5;
	}
	else
	{
	points_length2 = 2;
	}
	if (fil2->height * 2 > fil2->width)
	{
	points_height2 = 5;
	}
	else
	{
	points_height2 = 2;
	}
	}
	else
	{
	points_height2 = 5;

	if (fil2->length * 2 > fil2->height)
	{
	points_length2 = 5;
	}
	else
	{
	points_length2 = 2;
	}
	if (fil2->width * 2 > fil2->height)
	{
	points_width2 = 5;
	}
	else
	{
	points_width2 = 2;
	}

	}
	*/


	double *wi00 = (double *)malloc(points_length2 * sizeof(double));
	double *wi11 = (double *)malloc(points_width2 * sizeof(double));
	double *wi22 = (double *)malloc(points_height2 * sizeof(double));

	double *xi00 = (double *)malloc(points_length2 * sizeof(double));
	double *xi11 = (double *)malloc(points_width2 * sizeof(double));
	double *xi22 = (double *)malloc(points_height2 * sizeof(double));

	double *length_vec2 = (double *)malloc(points_length2 * 3 * sizeof(double));
	double *width_vec2 = (double *)malloc(points_width2 * 3 * sizeof(double));
	double *height_vec2 = (double *)malloc(points_height2 * 3 * sizeof(double));


	if (points4 == 1)
	{
		wi00[0] = 1.0000000000000000e+000;
	}
	else if (points4 == 2)
	{
		wi00[0] = 4.9999999999999989e-001;
		wi00[1] = 4.9999999999999989e-001;
	}
	else if (points4 == 3)
	{
		wi00[0] = 2.7777777777777768e-001;
		wi00[1] = 4.4444444444444425e-001;
		wi00[2] = 2.7777777777777779e-001;

	}
	else if (points4 == 4)
	{
		wi00[0] = 1.7392742256872706e-001;
		wi00[1] = 3.2607257743127316e-001;
		wi00[2] = 3.2607257743127299e-001;
		wi00[3] = 1.7392742256872701e-001;
	}
	else if (points4 == 5)
	{
		wi00[0] = 1.1846344252809456e-001;
		wi00[1] = 2.3931433524968337e-001;
		wi00[2] = 2.8444444444444444e-001;
		wi00[3] = 2.3931433524968343e-001;
		wi00[4] = 1.1846344252809445e-001;


	}
	else if (points4 == 6)
	{
		wi00[0] = 8.5662246189585289e-002;
		wi00[1] = 1.8038078652406928e-001;
		wi00[2] = 2.3395696728634516e-001;
		wi00[3] = 2.3395696728634591e-001;
		wi00[4] = 1.8038078652406928e-001;
		wi00[5] = 8.5662246189585220e-002;
	}
	else if (points4 == 7)
	{
		wi00[0] = 6.4742483084434824e-002;
		wi00[1] = 1.3985269574463838e-001;
		wi00[2] = 1.9091502525255943e-001;
		wi00[3] = 2.0897959183673431e-001;
		wi00[4] = 1.9091502525255952e-001;
		wi00[5] = 1.3985269574463827e-001;
		wi00[6] = 6.4742483084434907e-002;
	}


	if (points5 == 1)
	{
		wi11[0] = 1.0000000000000000e+000;
	}
	else if (points5 == 2)
	{
		wi11[0] = 4.9999999999999989e-001;
		wi11[1] = 4.9999999999999989e-001;
	}
	else if (points5 == 3)
	{
		wi11[0] = 2.7777777777777768e-001;
		wi11[1] = 4.4444444444444425e-001;
		wi11[2] = 2.7777777777777779e-001;
	}
	else if (points5 == 4)
	{
		wi11[0] = 1.7392742256872706e-001;
		wi11[1] = 3.2607257743127316e-001;
		wi11[2] = 3.2607257743127299e-001;
		wi11[3] = 1.7392742256872701e-001;
	}
	else if (points5 == 5)
	{
		wi11[0] = 1.1846344252809456e-001;
		wi11[1] = 2.3931433524968337e-001;
		wi11[2] = 2.8444444444444444e-001;
		wi11[3] = 2.3931433524968343e-001;
		wi11[4] = 1.1846344252809445e-001;
	}
	else if (points5 == 6)
	{
		wi11[0] = 8.5662246189585289e-002;
		wi11[1] = 1.8038078652406928e-001;
		wi11[2] = 2.3395696728634516e-001;
		wi11[3] = 2.3395696728634591e-001;
		wi11[4] = 1.8038078652406928e-001;
		wi11[5] = 8.5662246189585220e-002;
	}
	else if (points5 == 7)
	{
		wi11[0] = 6.4742483084434824e-002;
		wi11[1] = 1.3985269574463838e-001;
		wi11[2] = 1.9091502525255943e-001;
		wi11[3] = 2.0897959183673431e-001;
		wi11[4] = 1.9091502525255952e-001;
		wi11[5] = 1.3985269574463827e-001;
		wi11[6] = 6.4742483084434907e-002;
	}

	if (points6 == 1)
	{
		wi22[0] = 1.0000000000000000e+000;
	}
	else if (points6 == 2)
	{
		wi22[0] = 4.9999999999999989e-001;
		wi22[1] = 4.9999999999999989e-001;
	}
	else if (points6 == 3)
	{
		wi22[0] = 2.7777777777777768e-001;
		wi22[1] = 4.4444444444444425e-001;
		wi22[2] = 2.7777777777777779e-001;

	}
	else if (points6 == 4)
	{
		wi22[0] = 1.7392742256872706e-001;
		wi22[1] = 3.2607257743127316e-001;
		wi22[2] = 3.2607257743127299e-001;
		wi22[3] = 1.7392742256872701e-001;
	}
	else if (points6 == 5)
	{
		wi22[0] = 1.1846344252809456e-001;
		wi22[1] = 2.3931433524968337e-001;
		wi22[2] = 2.8444444444444444e-001;
		wi22[3] = 2.3931433524968343e-001;
		wi22[4] = 1.1846344252809445e-001;
	}
	else if (points6 == 6)
	{
		wi22[0] = 8.5662246189585289e-002;
		wi22[1] = 1.8038078652406928e-001;
		wi22[2] = 2.3395696728634516e-001;
		wi22[3] = 2.3395696728634591e-001;
		wi22[4] = 1.8038078652406928e-001;
		wi22[5] = 8.5662246189585220e-002;
	}
	else if (points6 == 7)
	{

		wi22[0] = 6.4742483084434824e-002;
		wi22[1] = 1.3985269574463838e-001;
		wi22[2] = 1.9091502525255943e-001;
		wi22[3] = 2.0897959183673431e-001;
		wi22[4] = 1.9091502525255952e-001;
		wi22[5] = 1.3985269574463827e-001;
		wi22[6] = 6.4742483084434907e-002;

	}



	if (points4 == 1)
	{
		xi00[0] = 5.0000000000000000e-001;
	}
	else if (points4 == 2)
	{
		xi00[0] = 2.1132486540518708e-001;
		xi00[1] = 7.8867513459481287e-001;
	}
	else if (points4 == 3)
	{
		xi00[0] = 1.1270166537925824e-001;
		xi00[1] = 5.0000000000000000e-001;
		xi00[2] = 8.8729833462074170e-001;

	}
	else if (points4 == 4)
	{
		xi00[0] = 6.9431844202973714e-002;
		xi00[1] = 3.3000947820757182e-001;
		xi00[2] = 6.6999052179242813e-001;
		xi00[3] = 9.3056815579702634e-001;
	}
	else if (points4 == 5)
	{
		xi00[0] = 4.6910077030668074e-002;
		xi00[1] = 2.3076534494715845e-001;
		xi00[2] = 5.0000000000000000e-001;
		xi00[3] = 7.6923465505284150e-001;
		xi00[4] = 9.5308992296933193e-001;
	}
	else if (points4 == 6)
	{
		xi00[0] = 3.3765242898423975e-002;
		xi00[1] = 1.6939530676686776e-001;
		xi00[2] = 3.8069040695840151e-001;
		xi00[3] = 6.1930959304159849e-001;
		xi00[4] = 8.3060469323313224e-001;
		xi00[5] = 9.6623475710157614e-001;
	}
	else if (points4 == 7)
	{

		xi00[0] = 2.5446043828620812e-002;
		xi00[1] = 1.2923440720030277e-001;
		xi00[2] = 2.9707742431130146e-001;
		xi00[3] = 4.9999999999999994e-001;
		xi00[4] = 7.0292257568869854e-001;
		xi00[5] = 8.7076559279969734e-001;
		xi00[6] = 9.7455395617137919e-001;
	}

	if (points5 == 1)
	{
		xi11[0] = 5.0000000000000000e-001;
	}
	else if (points5 == 2)
	{
		xi11[0] = 2.1132486540518708e-001;
		xi11[1] = 7.8867513459481287e-001;
	}
	else if (points5 == 3)
	{
		xi11[0] = 1.1270166537925824e-001;
		xi11[1] = 5.0000000000000000e-001;
		xi11[2] = 8.8729833462074170e-001;
	}
	else if (points5 == 4)
	{
		xi11[0] = 6.9431844202973714e-002;
		xi11[1] = 3.3000947820757182e-001;
		xi11[2] = 6.6999052179242813e-001;
		xi11[3] = 9.3056815579702634e-001;
	}
	else if (points5 == 5)
	{
		xi11[0] = 4.6910077030668074e-002;
		xi11[1] = 2.3076534494715845e-001;
		xi11[2] = 5.0000000000000000e-001;
		xi11[3] = 7.6923465505284150e-001;
		xi11[4] = 9.5308992296933193e-001;
	}
	else if (points5 == 6)
	{
		xi11[0] = 3.3765242898423975e-002;
		xi11[1] = 1.6939530676686776e-001;
		xi11[2] = 3.8069040695840151e-001;
		xi11[3] = 6.1930959304159849e-001;
		xi11[4] = 8.3060469323313224e-001;
		xi11[5] = 9.6623475710157614e-001;
	}
	else if (points5 == 7)
	{

		xi11[0] = 2.5446043828620812e-002;
		xi11[1] = 1.2923440720030277e-001;
		xi11[2] = 2.9707742431130146e-001;
		xi11[3] = 4.9999999999999994e-001;
		xi11[4] = 7.0292257568869854e-001;
		xi11[5] = 8.7076559279969734e-001;
		xi11[6] = 9.7455395617137919e-001;
	}




	if (points6 == 1)
	{
		xi22[0] = 5.0000000000000000e-001;
	}
	else if (points6 == 2)
	{
		xi22[0] = 2.1132486540518708e-001;
		xi22[1] = 7.8867513459481287e-001;
	}
	else if (points6 == 3)
	{
		xi22[0] = 1.1270166537925824e-001;
		xi22[1] = 5.0000000000000000e-001;
		xi22[2] = 8.8729833462074170e-001;
	}
	else if (points6 == 4)
	{
		xi22[0] = 6.9431844202973714e-002;
		xi22[1] = 3.3000947820757182e-001;
		xi22[2] = 6.6999052179242813e-001;
		xi22[3] = 9.3056815579702634e-001;
	}
	else if (points6 == 5)
	{
		xi22[0] = 4.6910077030668074e-002;
		xi22[1] = 2.3076534494715845e-001;
		xi22[2] = 5.0000000000000000e-001;
		xi22[3] = 7.6923465505284150e-001;
		xi22[4] = 9.5308992296933193e-001;
	}
	else if (points6 == 6)
	{
		xi22[0] = 3.3765242898423975e-002;
		xi22[1] = 1.6939530676686776e-001;
		xi22[2] = 3.8069040695840151e-001;
		xi22[3] = 6.1930959304159849e-001;
		xi22[4] = 8.3060469323313224e-001;
		xi22[5] = 9.6623475710157614e-001;
	}
	else if (points6 == 7)
	{
		xi22[0] = 2.5446043828620812e-002;
		xi22[1] = 1.2923440720030277e-001;
		xi22[2] = 2.9707742431130146e-001;
		xi22[3] = 4.9999999999999994e-001;
		xi22[4] = 7.0292257568869854e-001;
		xi22[5] = 8.7076559279969734e-001;
		xi22[6] = 9.7455395617137919e-001;
	}


	for (int i = 0; i < points4; i++)
	{
		wi00[i] *= 2;

	}

	for (int i = 0; i < points5; i++)
	{
		wi11[i] *= 2;

	}

	for (int i = 0; i < points6; i++)
	{
		wi22[i] *= 2;

	}

	for (int i = 0; i < points4; i++)
	{
		xi00[i] = xi00[i] * 2 - 1;

	}

	for (int i = 0; i < points5; i++)
	{
		xi11[i] = xi11[i] * 2 - 1;

	}


	for (int i = 0; i < points6; i++)
	{
		xi22[i] = xi22[i] * 2 - 1;

	}
	/*
	if (points_length2 == 5)
	{
	wi00[0] = 0.5688888888888889;
	wi00[1] = 0.4786286704993665;
	wi00[2] = 0.4786286704993665;
	wi00[3] = 0.2369268850561891;
	wi00[4] = 0.2369268850561891;

	xi00[0] = 0;
	xi00[1] = -0.5384693101056831;
	xi00[2] = 0.5384693101056831;
	xi00[3] = -0.9061798459386640;
	xi00[4] = 0.9061798459386640;
	}
	else
	{
	wi00[0] = 1;
	wi00[1] = 1;

	xi00[0] = -0.5773502691896257;
	xi00[1] = 0.5773502691896257;
	}

	*/
	for (int i = 0; i < points_length2; i++)
	{
		length_vec2[i] = (((fil2->x[1] - fil2->x[0]) / 2)*xi00[i])*(fil2->lenvect[0] / (sqrt(pow(fil2->lenvect[0], 2) + pow(fil2->lenvect[1], 2) + pow(fil2->lenvect[2], 2))));
		length_vec2[i + points_length2] = (((fil2->y[1] - fil2->y[0]) / 2)*xi00[i])*(fil2->lenvect[1] / (sqrt(pow(fil2->lenvect[0], 2) + pow(fil2->lenvect[1], 2) + pow(fil2->lenvect[2], 2))));
		length_vec2[i + points_length2 * 2] = (((fil2->z[1] - fil2->z[0]) / 2)*xi00[i])*(fil2->lenvect[2] / (sqrt(pow(fil2->lenvect[0], 2) + pow(fil2->lenvect[1], 2) + pow(fil2->lenvect[2], 2))));
	}

	/*


	if (points_width2 == 5)
	{
	wi11[0] = 0.5688888888888889;
	wi11[1] = 0.4786286704993665;
	wi11[2] = 0.4786286704993665;
	wi11[3] = 0.2369268850561891;
	wi11[4] = 0.2369268850561891;

	xi11[0] = 0;
	xi11[1] = -0.5384693101056831;
	xi11[2] = 0.5384693101056831;
	xi11[3] = -0.9061798459386640;
	xi11[4] = 0.9061798459386640;
	}
	else
	{
	wi11[0] = 1;
	wi11[1] = 1;

	xi11[0] = -0.5773502691896257;
	xi11[1] = 0.5773502691896257;
	}
	*/
	for (int i = 0; i < points_width2; i++)
	{
		width_vec2[i] = (((fil2->width) / 2)*xi11[i])*width2[0];
		width_vec2[i + points_width2] = (((fil2->width) / 2)*xi11[i])*width2[1];
		width_vec2[i + points_width2 * 2] = (((fil2->width) / 2)*xi11[i])*width2[2];
	}

	/*
	if (points_height2 == 5)
	{
	wi22[0] = 0.5688888888888889;
	wi22[1] = 0.4786286704993665;
	wi22[2] = 0.4786286704993665;
	wi22[3] = 0.2369268850561891;
	wi22[4] = 0.2369268850561891;

	xi22[0] = 0;
	xi22[1] = -0.5384693101056831;
	xi22[2] = 0.5384693101056831;
	xi22[3] = -0.9061798459386640;
	xi22[4] = 0.9061798459386640;
	}
	else
	{
	wi22[0] = 1;
	wi22[1] = 1;

	xi22[0] = -0.5773502691896257;
	xi22[1] = 0.5773502691896257;
	}
	*/
	for (int i = 0; i < points_height2; i++)
	{
		height_vec2[i] = (((fil2->height) / 2)*xi22[i])*height2[0];
		height_vec2[i + points_height2] = (((fil2->height) / 2)*xi22[i])*height2[1];
		height_vec2[i + points_height2 * 2] = (((fil2->height) / 2)*xi22[i])*height2[2];

	}



	//height1[0] = fil1->lenvect[1] * width1[2] - fil1->lenvect[2] * width1[1];
	//height1[1] = fil1->lenvect[2] * width1[0] - fil1->lenvect[0] * width1[2];
	//height1[2] = fil1->lenvect[0] * width1[1] - fil1->lenvect[1] * width1[0];

	//double magg = sqrt(pow(height1[0],2) + pow(height1[1],2) + pow(height1[2],2));

	//height1[0] = height1[0]/magg;
	//height1[1] = height1[1]/magg;
	//height1[2] = height1[2]/magg;

	//height2[0] = fil2->lenvect[1] * width2[2] - fil2->lenvect[2] * width2[1];
	//height2[1] = fil2->lenvect[2] * width2[0] - fil2->lenvect[0] * width2[2];
	//height2[2] = fil2->lenvect[0] * width2[1] - fil2->lenvect[1] * width2[0];

	//magg = sqrt(pow(height2[0], 2) + pow(height2[1], 2) + pow(height2[2], 2));

	//height2[0] = height2[0] / magg;
	//height2[1] = height2[1] / magg;
	//height2[2] = height2[2] / magg;
	/*

	for (int i = 0; i < points; i++)
	{
	*/
	/*
	length_vec1[i] = (((fil1->x[1] - fil1->x[0]) / 2)*xi[i] + (fil1->x[1] + fil1->x[0])/2)*(fil1->lenvect[0]);
	length_vec1[i + 5] = (((fil1->y[1] - fil1->y[0]) / 2)*xi[i] + (fil1->y[1] + fil1->y[0]) / 2);
	length_vec1[i + 10] = (((fil1->z[1] - fil1->z[0]) / 2)*xi[i] + (fil1->z[1] + fil1->z[0]) / 2);

	length_vec2[i] = (((fil2->x[1] - fil2->x[0]) / 2)*xi[i] + (fil2->x[1] + fil2->x[0]) / 2);
	length_vec2[i + 5] = (((fil2->y[1] - fil2->y[0]) / 2)*xi[i] + (fil2->y[1] + fil2->y[0]) / 2);
	length_vec1[i + 10] = (((fil2->z[1] - fil2->z[0]) / 2)*xi[i] + (fil2->z[1] + fil2->z[0]) / 2);
	*/
	/*
	length_vec1[i] = (((fil1->x[1] - fil1->x[0]) / 2)*xi[i])*(fil1->lenvect[0]/(sqrt(pow(fil1->lenvect[0], 2) + pow(fil1->lenvect[1], 2) + pow(fil1->lenvect[2], 2))));
	length_vec1[i + points] = (((fil1->y[1] - fil1->y[0]) / 2)*xi[i])*(fil1->lenvect[1] / (sqrt(pow(fil1->lenvect[0], 2) + pow(fil1->lenvect[1], 2) + pow(fil1->lenvect[2], 2))));
	length_vec1[i + points*2] = (((fil1->z[1] - fil1->z[0]) / 2)*xi[i])*(fil1->lenvect[2] / (sqrt(pow(fil1->lenvect[0], 2) + pow(fil1->lenvect[1], 2) + pow(fil1->lenvect[2], 2))));

	length_vec2[i] = (((fil2->x[1] - fil2->x[0]) / 2)*xi[i])*(fil2->lenvect[0] / (sqrt(pow(fil2->lenvect[0], 2) + pow(fil2->lenvect[1], 2) + pow(fil2->lenvect[2], 2))));
	length_vec2[i + points] = (((fil2->y[1] - fil2->y[0]) / 2)*xi[i])*(fil2->lenvect[1] / (sqrt(pow(fil2->lenvect[0], 2) + pow(fil2->lenvect[1], 2) + pow(fil2->lenvect[2], 2))));
	length_vec2[i + points*2] = (((fil2->z[1] - fil2->z[0]) / 2)*xi[i])*(fil2->lenvect[2] / (sqrt(pow(fil2->lenvect[0], 2) + pow(fil2->lenvect[1], 2) + pow(fil2->lenvect[2], 2))));
	}
	*/
	/*
	for (int i = 0; i < edge; i++)
	{
	*/
	/*
	width_vec1[i] = (((fil1->width) / 2)*xi2[i] + fil1->x[0])*width1[0];
	width_vec1[i+2] = (((fil1->width) / 2)*xi2[i] + fil1->y[0])*width1[1];
	width_vec1[i+4] = (((fil1->width) / 2)*xi2[i] + fil1->z[0])*width1[2];

	width_vec2[i] = (((fil2->width) / 2)*xi2[i] + fil2->x[0])*width2[0];
	width_vec2[i + 2] = (((fil2->width) / 2)*xi2[i] + fil2->y[0])*width2[1];
	width_vec2[i + 4] = (((fil2->width) / 2)*xi2[i] + fil2->z[0])*width2[2];

	height_vec1[i] = (((fil1->height) / 2)*xi2[i] + fil1->x[0])*height1[0];
	height_vec1[i+2] = (((fil1->height) / 2)*xi2[i] + fil1->y[0])*height1[1];
	height_vec1[i+4] = (((fil1->height) / 2)*xi2[i] + fil1->z[0])*height1[2];


	height_vec2[i] = (((fil2->height) / 2)*xi2[i] + fil2->x[0])*height2[0];
	height_vec2[i + 2] = (((fil2->height) / 2)*xi2[i] + fil2->y[0])*height2[1];
	height_vec2[i + 4] = (((fil2->height) / 2)*xi2[i] + fil2->z[0])*height2[2];
	*/
	/*
	width_vec1[i] = (((fil1->width) / 2)*xi2[i])*width1[0];
	width_vec1[i + edge] = (((fil1->width) / 2)*xi2[i])*width1[1];
	width_vec1[i + edge*2] = (((fil1->width) / 2)*xi2[i])*width1[2];

	width_vec2[i] = (((fil2->width) / 2)*xi2[i])*width2[0];
	width_vec2[i + edge] = (((fil2->width) / 2)*xi2[i])*width2[1];
	width_vec2[i + edge*2] = (((fil2->width) / 2)*xi2[i])*width2[2];

	height_vec1[i] = (((fil1->height) / 2)*xi2[i])*height1[0];
	height_vec1[i + edge] = (((fil1->height) / 2)*xi2[i])*height1[1];
	height_vec1[i + edge*2] = (((fil1->height) / 2)*xi2[i])*height1[2];


	height_vec2[i] = (((fil2->height) / 2)*xi2[i])*height2[0];
	height_vec2[i + edge] = (((fil2->height) / 2)*xi2[i])*height2[1];
	height_vec2[i + edge*2] = (((fil2->height) / 2)*xi2[i])*height2[2];

	}
	*/


	double nnn, mmm, lll, kkk, jjj, iii = 0;
	double x11, y11, z11, x22, y22, z22;

	for (int i = 0; i < points_length; i++)
	{
		jjj = 0;
		for (int j = 0; j < points_width; j++)
		{
			kkk = 0;
			for (int k = 0; k < points_height; k++)
			{
				lll = 0;
				x11 = height_vec1[k] + width_vec1[j] + length_vec1[i] + (fil1->x[0] + fil1->x[1]) / 2;
				y11 = height_vec1[k + points_height] + width_vec1[j + points_width] + length_vec1[i + points_length] + (fil1->y[0] + fil1->y[1]) / 2;
				z11 = height_vec1[k + points_height * 2] + width_vec1[j + points_width * 2] + length_vec1[i + points_length * 2] + (fil1->z[0] + fil1->z[1]) / 2;
				for (int l = 0; l < points_length2; l++)
				{
					mmm = 0;
					for (int m = 0; m < points_width2; m++)
					{
						nnn = 0;
						for (int n = 0; n < points_height2; n++)
						{
							x22 = height_vec2[n] + width_vec2[m] + length_vec2[l] + (fil2->x[0] + fil2->x[1]) / 2;
							y22 = height_vec2[n + points_height2] + width_vec2[m + points_width2] + length_vec2[l + points_length2] + (fil2->y[0] + fil2->y[1]) / 2;
							z22 = height_vec2[n + points_height2 * 2] + width_vec2[m + points_width2 * 2] + length_vec2[l + points_length2 * 2] + (fil2->z[0] + fil2->z[1]) / 2;
							nnn += wi22[n] * (1 / sqrt(pow(x11 - x22, 2) + pow(y11 - y22, 2) + pow(z11 - z22, 2)));
						}
						mmm += wi11[m] * nnn;
					}
					lll += wi00[l] * mmm;
				}
				kkk += wi2[k] * lll;
			}
			jjj += wi1[j] * kkk;
		}
		iii += wi[i] * jjj;
	}

	//return  (1) * a_b * iii;

	free(length_vec1);
	free(length_vec2);

	free(width_vec1);
	free(width_vec2);

	free(height_vec1);
	free(height_vec2);

	free(wi);
	free(wi1);
	free(wi2);
	free(wi00);
	free(wi11);
	free(wi22);

	free(xi);
	free(xi1);
	free(xi2);
	free(xi00);
	free(xi11);
	free(xi22);
	return  iii / 64;
}


double ind(FILAMENT *fil11, FILAMENT *fil22, double sum)
{
	double dist = sqrt(pow((fil11->x[0] + fil11->x[1]) / 2 - (fil22->x[0] + fil22->x[1]) / 2, 2)
		+ pow((fil11->y[0] + fil11->y[1]) / 2 - (fil22->y[0] + fil22->y[1]) / 2, 2)
		+ pow((fil11->z[0] + fil11->z[1]) / 2 - (fil22->z[0] + fil22->z[1]) / 2, 2));
	double length;
	if (fil11->max_dim > fil22->max_dim)
	{
		length = fil11->max_dim;
	}
	else
	{
		length = fil22->max_dim;
	}



	if (41 * length < dist)
	{
		return	mutual_quadrature(fil11, fil22, 1, 1, 1, 1, 1, 1) - sum;
	}
	/*else if ((fil11->x_min - 4*length) > fil22->x_max || (fil11->x_max + 4*length) <  fil22->x_min
	|| (fil11->y_min - 4*length) > fil22->y_max || (fil11->y_max + 4*length) <  fil22->y_min
	|| (fil11->z_min - 4*length) > fil22->z_max || (fil11->z_max + 4*length) <  fil22->z_min || true)*/
	else if (5 * length < dist)
	{


		if (fil11->max_dim == fil11->length)
		{
			if (fil11->max_dim == fil11->length)
			{
				return	mutual_quadrature(fil11, fil22, 1, 2, 2, 1, 2, 2) - sum;
			}
			else if (fil11->max_dim == fil11->width)
			{
				return	mutual_quadrature(fil11, fil22, 1, 2, 2, 2, 1, 2) - sum;
			}
			else if (fil11->max_dim == fil11->height)
			{
				return	mutual_quadrature(fil11, fil22, 1, 2, 2, 2, 2, 1) - sum;
			}
			else
			{
				return	mutual_quadrature(fil11, fil22, 1, 2, 2, 2, 2, 2) - sum;
			}
		}
		else if (fil11->max_dim == fil11->width)
		{
			if (fil11->max_dim == fil11->length)
			{
				return	mutual_quadrature(fil11, fil22, 2, 1, 2, 1, 2, 2) - sum;
			}
			else if (fil11->max_dim == fil11->width)
			{
				return	mutual_quadrature(fil11, fil22, 2, 1, 2, 2, 1, 2) - sum;
			}
			else if (fil11->max_dim == fil11->height)
			{
				return	mutual_quadrature(fil11, fil22, 2, 1, 2, 2, 2, 1) - sum;
			}
			else
			{
				return	mutual_quadrature(fil11, fil22, 2, 1, 2, 2, 2, 2) - sum;
			}
		}
		else if (fil11->max_dim == fil11->height)
		{
			if (fil11->max_dim == fil11->length)
			{
				return	mutual_quadrature(fil11, fil22, 2, 2, 1, 1, 2, 2) - sum;
			}
			else if (fil11->max_dim == fil11->width)
			{
				return	mutual_quadrature(fil11, fil22, 2, 2, 1, 2, 1, 2) - sum;
			}
			else if (fil11->max_dim == fil11->height)
			{
				return	mutual_quadrature(fil11, fil22, 2, 2, 1, 2, 2, 1) - sum;
			}
			else
			{
				return	mutual_quadrature(fil11, fil22, 2, 2, 1, 2, 2, 2) - sum;
			}
		}
		else
		{
			return	mutual_quadrature(fil11, fil22, 2, 2, 2, 2, 2, 2) - sum;
		}
	}

	/*else if ((fil11->x_min - 3 * length) > fil22->x_max || (fil11->x_max + 3 * length) <  fil22->x_min
	|| (fil11->y_min - 3 * length) > fil22->y_max || (fil11->y_max + 3 * length) <  fil22->y_min
	|| (fil11->z_min - 3 * length) > fil22->z_max || (fil11->z_max + 3 * length) <  fil22->z_min)*/
	else if (4 * length < dist)
	{
		return	mutual_quadrature(fil11, fil22, 2, 2, 2, 2, 2, 2) - sum;
	}

	else
	{

		//first test seperation 
		if ((fil11->x_min - 0.15*length) < fil22->x_max && (fil11->x_max + 0.15*length) > fil22->x_min
			&& (fil11->y_min - 0.15*length) < fil22->y_max && (fil11->y_max + 0.15*length) > fil22->y_min
			&& (fil11->z_min - 0.15*length) < fil22->z_max && (fil11->z_max + 0.15*length) > fil22->z_min)
		{




			if (fil11->min_dim == fil11->length)
			{
				if (fil11->min_dim == fil11->width && fil11->min_dim == fil11->height)
				{
					return aanalitic(fil11, fil22, 4, 4, 4) - sum;
				}
				else if (fil11->min_dim == fil11->width)
				{
					if (fil11->min_dim * 10 > fil11->height)
					{
						return aanalitic(fil11, fil22, 4, 4, 8) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 4, 4, 18) - sum;
					}
				}
				else if (fil11->min_dim == fil11->height)
				{
					if (fil11->min_dim * 10 > fil11->width)
					{
						return aanalitic(fil11, fil22, 4, 8, 4) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 4, 18, 4) - sum;
					}
				}
				else if (fil11->min_dim * 10 > fil11->width)
				{
					if (fil11->min_dim * 10 > fil11->height)
					{
						return aanalitic(fil11, fil22, 4, 8, 8) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 4, 8, 18) - sum;
					}
				}
				else
				{
					if (fil11->min_dim * 10 < fil11->height)
					{
						return aanalitic(fil11, fil22, 4, 18, 18) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 4, 18, 8) - sum;
					}
				}


			}


			else if (fil11->min_dim == fil11->height)
			{
				if (fil11->min_dim == fil11->width)
				{
					if (fil11->min_dim * 10 > fil11->length)
					{
						return aanalitic(fil11, fil22, 8, 4, 4) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 18, 4, 4) - sum;
					}
				}
				else if (fil11->min_dim * 10 > fil11->width)
				{
					if (fil11->min_dim * 10 > fil11->length)
					{
						return aanalitic(fil11, fil22, 8, 8, 4) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 18, 8, 4) - sum;
					}
				}
				else
				{
					if (fil11->min_dim * 10 < fil11->length)
					{
						return aanalitic(fil11, fil22, 18, 18, 4) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 8, 18, 4) - sum;
					}
				}
			}


			else if (fil11->min_dim == fil11->width)
			{



				if (fil11->min_dim * 10 > fil11->height)
				{
					if (fil11->min_dim * 10 > fil11->length)
					{
						return aanalitic(fil11, fil22, 8, 4, 8) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 18, 4, 8) - sum;
					}
				}
				else
				{
					if (fil11->min_dim * 10 < fil11->length)
					{
						return aanalitic(fil11, fil22, 18, 4, 18) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 8, 4, 18) - sum;
					}
				}
			}


		}


		if ((fil11->x_min - 0.3*length) < fil22->x_max && (fil11->x_max + 0.3*length) > fil22->x_min
			&& (fil11->y_min - 0.3*length) < fil22->y_max && (fil11->y_max + 0.3*length) > fil22->y_min
			&& (fil11->z_min - 0.3*length) < fil22->z_max && (fil11->z_max + 0.3*length) > fil22->z_min)
		{
			if (fil11->min_dim == fil11->length)
			{

				return aanalitic(fil11, fil22, 3, 5, 5) - sum;

			}


			else if (fil11->min_dim == fil11->height)
			{
				return aanalitic(fil11, fil22, 5, 5, 3) - sum;
			}


			else if (fil11->min_dim == fil11->width)
			{


				return aanalitic(fil11, fil22, 5, 3, 5) - sum;

			}
		}



		/*else if ((fil11->x_min - 0.3 * length) > fil22->x_max || (fil11->x_max + 0.3 * length) < fil22->x_min
		|| (fil11->y_min - 0.3 * length) > fil22->y_max || (fil11->y_max + 0.3 * length) < fil22->y_min
		|| (fil11->z_min - 0.3 * length) > fil22->z_max || (fil11->z_max + 0.3 * length) < fil22->z_min)*/
		else
		{
			int x1, y1, z1, x2, y2, z2;

			if (fil11->length == fil11->min_dim)
			{
				x1 = 3;
			}
			else
			{
				x1 = 6;
			}

			if (fil11->width == fil11->min_dim)
			{
				y1 = 3;
			}
			else
			{
				y1 = 6;
			}
			if (fil11->height == fil11->min_dim)
			{
				z1 = 3;
			}
			else
			{
				z1 = 6;
			}


			if (fil11->length > fil11->min_dim * 10)
			{
				x1 = 7;
			}
			if (fil11->width > fil11->min_dim * 10)
			{
				y1 = 7;
			}
			if (fil11->height > fil11->min_dim * 10)
			{
				z1 = 7;
			}



			if (fil22->length == fil22->min_dim)
			{
				x2 = 3;
			}
			else
			{
				x2 = 6;
			}

			if (fil22->width == fil22->min_dim)
			{
				y2 = 3;
			}
			else
			{
				y2 = 6;
			}
			if (fil22->height == fil22->min_dim)
			{
				z2 = 3;
			}
			else
			{
				z2 = 6;
			}


			if (fil22->length > fil22->min_dim * 10)
			{
				x2 = 7;
			}
			if (fil22->width > fil22->min_dim * 10)
			{
				y2 = 7;
			}
			if (fil22->height > fil22->min_dim * 10)
			{
				z2 = 7;
			}


			return	mutual_quadrature(fil11, fil22, x1, y1, z1, x2, y2, z2) - sum;
		}

	}
	return 0;
}













/*Determines the number of integration points to be used and 
whether Aanalitic or mutual_quadrature is required.*/
double induct(FILAMENT *fil11, FILAMENT *fil22, double sum)
{

	double dist = sqrt(pow((fil11->x[0] + fil11->x[1]) / 2 - (fil22->x[0] + fil22->x[1]) / 2, 2)
		+ pow((fil11->y[0] + fil11->y[1]) / 2 - (fil22->y[0] + fil22->y[1]) / 2, 2)
		+ pow((fil11->z[0] + fil11->z[1]) / 2 - (fil22->z[0] + fil22->z[1]) / 2, 2));
	double length;
	if (fil11->max_dim > fil22->max_dim)
	{
		length = fil11->max_dim;
	}
	else
	{
		length = fil22->max_dim;
	}

	if (450 * length < dist)
	{
		return	mutual_quadrature(fil11, fil22, 1, 1, 1, 1, 1, 1) - sum;
	}
	else if (12 * length < dist)
	{
		return	mutual_quadrature(fil11, fil22, 2, 2, 2, 2, 2, 2) - sum;
	}
	if (2.5*length < dist)
	{
		if (fil11->max_dim == fil11->length)
		{
			if (fil11->max_dim == fil11->length)
			{
				return	mutual_quadrature(fil11, fil22, 3, 4, 4, 3, 4, 4) - sum;
			}
			else if (fil11->max_dim == fil11->width)
			{
				return	mutual_quadrature(fil11, fil22, 3, 4, 4, 4, 3, 4) - sum;
			}
			else if (fil11->max_dim == fil11->height)
			{
				return	mutual_quadrature(fil11, fil22, 3, 4, 4, 4, 4, 3) - sum;
			}
			else
			{
				return	mutual_quadrature(fil11, fil22, 3, 4, 4, 4, 4, 4) - sum;
			}
		}
		else if (fil11->max_dim == fil11->width)
		{
			if (fil11->max_dim == fil11->length)
			{
				return	mutual_quadrature(fil11, fil22, 4, 3, 4, 3, 4, 4) - sum;
			}
			else if (fil11->max_dim == fil11->width)
			{
				return	mutual_quadrature(fil11, fil22, 4, 3, 4, 4, 3, 4) - sum;
			}
			else if (fil11->max_dim == fil11->height)
			{
				return	mutual_quadrature(fil11, fil22, 4, 3, 4, 4, 4, 3) - sum;
			}
			else
			{
				return	mutual_quadrature(fil11, fil22, 4, 3, 4, 4, 4, 4) - sum;
			}
		}
		else if (fil11->max_dim == fil11->height)
		{
			if (fil11->max_dim == fil11->length)
			{
				return	mutual_quadrature(fil11, fil22, 4, 4, 3, 3, 4, 4) - sum;
			}
			else if (fil11->max_dim == fil11->width)
			{
				return	mutual_quadrature(fil11, fil22, 4, 4, 3, 4, 3, 4) - sum;
			}
			else if (fil11->max_dim == fil11->height)
			{
				return	mutual_quadrature(fil11, fil22, 4, 4, 3, 4, 4, 3) - sum;
			}
			else
			{
				return	mutual_quadrature(fil11, fil22, 4, 4, 3, 4, 4, 4) - sum;
			}
		}
		else
		{
			return	mutual_quadrature(fil11, fil22, 4, 4, 4, 4, 4, 4) - sum;
		}

	}
	else if (1.5 * length < dist)
	{
		if (fil11->min_dim == fil11->length)
		{
			if (fil22->min_dim == fil22->length)
			{
				return	mutual_quadrature(fil11, fil22, 4, 5, 5, 4, 5, 5) - sum;
			}
			else if (fil22->min_dim == fil22->width)
			{
				return	mutual_quadrature(fil11, fil22, 4, 5, 5, 5, 4, 5) - sum;
			}
			else if (fil22->min_dim == fil22->height)
			{
				return	mutual_quadrature(fil11, fil22, 4, 5, 5, 5, 5, 4) - sum;
			}
			else
			{
				return	mutual_quadrature(fil11, fil22, 4, 5, 5, 5, 5, 5) - sum;
			}
		}
		else if (fil11->min_dim == fil11->width)
		{
			if (fil22->min_dim == fil22->length)
			{
				return	mutual_quadrature(fil11, fil22, 5, 4, 5, 4, 5, 5) - sum;
			}
			else if (fil22->min_dim == fil22->width)
			{
				return	mutual_quadrature(fil11, fil22, 5, 4, 5, 5, 4, 5) - sum;
			}
			else if (fil22->min_dim == fil22->height)
			{
				return	mutual_quadrature(fil11, fil22, 5, 4, 5, 5, 5, 4) - sum;
			}
			else
			{
				return	mutual_quadrature(fil11, fil22, 5, 4, 5, 5, 5, 5) - sum;
			}
		}
		else if (fil11->min_dim == fil11->height)
		{
			if (fil22->min_dim == fil22->length)
			{
				return	mutual_quadrature(fil11, fil22, 5, 5, 4, 4, 5, 5) - sum;
			}
			else if (fil22->min_dim == fil22->width)
			{
				return	mutual_quadrature(fil11, fil22, 5, 5, 4, 5, 4, 5) - sum;
			}
			else if (fil22->min_dim == fil22->height)
			{
				return	mutual_quadrature(fil11, fil22, 5, 5, 4, 5, 5, 4) - sum;
			}
			else
			{
				return	mutual_quadrature(fil11, fil22, 5, 5, 4, 5, 5, 5) - sum;
			}
		}
		else
		{
			return	mutual_quadrature(fil11, fil22, 5, 5, 5, 5, 5, 5) - sum;
		}
	}
	else
	{

		//first test seperation 

		if ((fil11->x_min - 0.1*length) < fil22->x_max && (fil11->x_max + 0.1*length) >  fil22->x_min
			&& (fil11->y_min - 0.1*length) < fil22->y_max && (fil11->y_max + 0.1*length) >  fil22->y_min
			&& (fil11->z_min - 0.1*length) < fil22->z_max && (fil11->z_max + 0.1*length) >  fil22->z_min)
		{
			if (fil11->min_dim == fil11->length)
			{
				if (fil11->min_dim == fil11->width && fil11->min_dim == fil11->height)
				{
					return aanalitic(fil11, fil22, 7, 7, 7) - sum;
				}
				else if (fil11->min_dim == fil11->width)
				{
					if (fil11->min_dim * 10 > fil11->height)
					{
						return aanalitic(fil11, fil22, 7, 7, 13) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 7, 7, 31) - sum;
					}
				}
				else if (fil11->min_dim == fil11->height)
				{
					if (fil11->min_dim * 10 > fil11->width)
					{
						return aanalitic(fil11, fil22, 7, 13, 7) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 7, 31, 7) - sum;
					}
				}
				else if (fil11->min_dim * 10 > fil11->width)
				{
					if (fil11->min_dim * 10 > fil11->height)
					{
						return aanalitic(fil11, fil22, 7, 13, 13) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 7, 13, 31) - sum;
					}
				}
				else
				{
					if (fil11->min_dim * 10 < fil11->height)
					{
						return aanalitic(fil11, fil22, 7, 31, 31) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 7, 31, 13) - sum;
					}
				}


			}


			else if (fil11->min_dim == fil11->height)
			{
				if (fil11->min_dim == fil11->width)
				{
					if (fil11->min_dim * 10 > fil11->length)
					{
						return aanalitic(fil11, fil22, 13, 7, 7) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 31, 7, 7) - sum;
					}
				}
				else if (fil11->min_dim * 10 > fil11->width)
				{
					if (fil11->min_dim * 10 > fil11->length)
					{
						return aanalitic(fil11, fil22, 13, 13, 7) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 31, 13, 7) - sum;
					}
				}
				else
				{
					if (fil11->min_dim * 10 < fil11->length)
					{
						return aanalitic(fil11, fil22, 31, 31, 7) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 13, 31, 7) - sum;
					}
				}
			}


			else if (fil11->min_dim == fil11->width)
			{



				if (fil11->min_dim * 10 > fil11->height)
				{
					if (fil11->min_dim * 10 > fil11->length)
					{
						return aanalitic(fil11, fil22, 13, 7, 13) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 31, 7, 13) - sum;
					}
				}
				else
				{
					if (fil11->min_dim * 10 < fil11->length)
					{
						return aanalitic(fil11, fil22, 31, 7, 31) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 13, 7, 31) - sum;
					}
				}
			}
		}
		else if ((fil11->x_min - 0.5*length) < fil22->x_max && (fil11->x_max + 0.5*length) > fil22->x_min
			&& (fil11->y_min - 0.5*length) < fil22->y_max && (fil11->y_max + 0.5*length) > fil22->y_min
			&& (fil11->z_min - 0.5*length) < fil22->z_max && (fil11->z_max + 0.5*length) > fil22->z_min)
		{




			if (fil11->min_dim == fil11->length)
			{
				if (fil11->min_dim == fil11->width && fil11->min_dim == fil11->height)
				{
					return aanalitic(fil11, fil22, 5, 5, 5) - sum;
				}
				else if (fil11->min_dim == fil11->width)
				{
					if (fil11->min_dim * 10 > fil11->height)
					{
						return aanalitic(fil11, fil22, 5, 5, 8) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 5, 5, 10) - sum;
					}
				}
				else if (fil11->min_dim == fil11->height)
				{
					if (fil11->min_dim * 10 > fil11->width)
					{
						return aanalitic(fil11, fil22, 5, 8, 5) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 5, 10, 5) - sum;
					}
				}
				else if (fil11->min_dim * 10 > fil11->width)
				{
					if (fil11->min_dim * 10 > fil11->height)
					{
						return aanalitic(fil11, fil22, 5, 8, 8) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 5, 8, 10) - sum;
					}
				}
				else
				{
					if (fil11->min_dim * 10 < fil11->height)
					{
						return aanalitic(fil11, fil22, 5, 10, 10) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 5, 10, 8) - sum;
					}
				}


			}


			else if (fil11->min_dim == fil11->height)
			{
				if (fil11->min_dim == fil11->width)
				{
					if (fil11->min_dim * 10 > fil11->length)
					{
						return aanalitic(fil11, fil22, 8, 5, 5) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 10, 5, 5) - sum;
					}
				}
				else if (fil11->min_dim * 10 > fil11->width)
				{
					if (fil11->min_dim * 10 > fil11->length)
					{
						return aanalitic(fil11, fil22, 8, 8, 5) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 10, 8, 5) - sum;
					}
				}
				else
				{
					if (fil11->min_dim * 10 < fil11->length)
					{
						return aanalitic(fil11, fil22, 10, 10, 5) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 8, 10, 5) - sum;
					}
				}
			}


			else if (fil11->min_dim == fil11->width)
			{



				if (fil11->min_dim * 10 > fil11->height)
				{
					if (fil11->min_dim * 10 > fil11->length)
					{
						return aanalitic(fil11, fil22, 8, 5, 8) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 10, 5, 8) - sum;
					}
				}
				else
				{
					if (fil11->min_dim * 10 < fil11->length)
					{
						return aanalitic(fil11, fil22, 10, 5, 10) - sum;
					}
					else
					{
						return aanalitic(fil11, fil22, 8, 5, 10) - sum;
					}
				}
			}


		}
		else
		{
			if (fil11->min_dim == fil11->length)
			{
				if (fil22->min_dim == fil22->length)
				{
					return	mutual_quadrature(fil11, fil22, 4, 6, 6, 4, 6, 6) - sum;
				}
				else if (fil22->min_dim == fil22->width)
				{
					return	mutual_quadrature(fil11, fil22, 4, 6, 6, 6, 4, 6) - sum;
				}
				else if (fil22->min_dim == fil22->height)
				{
					return	mutual_quadrature(fil11, fil22, 4, 6, 6, 6, 6, 4) - sum;
				}
				else
				{
					return	mutual_quadrature(fil11, fil22, 4, 6, 6, 6, 6, 6) - sum;
				}
			}
			else if (fil11->min_dim == fil11->width)
			{
				if (fil22->min_dim == fil22->length)
				{
					return	mutual_quadrature(fil11, fil22, 6, 4, 6, 4, 6, 6) - sum;
				}
				else if (fil22->min_dim == fil22->width)
				{
					return	mutual_quadrature(fil11, fil22, 6, 4, 6, 6, 4, 6) - sum;
				}
				else if (fil22->min_dim == fil22->height)
				{
					return	mutual_quadrature(fil11, fil22, 6, 4, 6, 6, 6, 4) - sum;
				}
				else
				{
					return	mutual_quadrature(fil11, fil22, 6, 4, 6, 6, 6, 6) - sum;
				}
			}
			else if (fil11->min_dim == fil11->height)
			{
				if (fil22->min_dim == fil22->length)
				{
					return	mutual_quadrature(fil11, fil22, 6, 6, 4, 4, 6, 6) - sum;
				}
				else if (fil22->min_dim == fil22->width)
				{
					return	mutual_quadrature(fil11, fil22, 6, 6, 4, 6, 4, 6) - sum;
				}
				else if (fil22->min_dim == fil22->height)
				{
					return	mutual_quadrature(fil11, fil22, 6, 6, 4, 6, 6, 4) - sum;
				}
				else
				{
					return	mutual_quadrature(fil11, fil22, 6, 6, 4, 6, 6, 6) - sum;
				}
			}
			else
			{
				return	mutual_quadrature(fil11, fil22, 6, 6, 6, 6, 6, 6) - sum;
			}
		}
	}
	return 0;
}



/*Computes the Frobenius norm of a vector*/
double frobeniusNormV3(double* matrix, int size1, int size2)
{
	double result = 0.0;
	for (int j = 0; j < size2; ++j)
	{
		//double value = *(matrix + (i*size2) + j);
		double value = matrix[j];
		result += value * value;
	}
	return sqrt(result);
}

/*
double frobeniusNormU3(double** matrix, int size1, int size2)
{
double result = 0.0;
for (int i = 0; i < size1; ++i)
{

//double value = **(matrix + (i*size2) + j);
double value = matrix[i][size2 - 1];
result += value * value;
}
return sqrt(result);
}
*/



/*Used to calculate a portion of Z ^ (k) for the recursive update and evaluation of convergence.*/
double multiplyy(int kk, int m, int n, double **U, double **V)
{


	double prodprod = 0;
	//sum
	for (int j = 0; j < kk; j++)//was k - 1
	{
		double prod1 = 0;
		double prod2 = 0;
		for (int i = 0; i < m; i++)//product1
		{
			//prod1 += abs(U(i, j)*U(i, kk));
			prod1 += abs(U[j][i] * U[kk][i]);
		}

		//product2
		for (int i = 0; i < n; i++)
		{
			//prod2 += abs(V(j, i)*V(kk, i));
			prod2 += abs(V[j][i] * V[kk][i]);
		}
		prodprod += 2 * prod1*prod2;
	}
	return prodprod;
}
/*
double  frobeniusUxV(double** matrixU, double* matrixV, int size1,int size2)
{
double **needed;
double **U = (double **)malloc(size1 * sizeof(double *));
for (int i = 0; i<size1; i++)
U[i] = (double *)malloc(size2 * sizeof(double));
for (int i = 0; i < size1; i++)
{

}
}
*/
//namespace ublas = boost::numeric::ublas;
//template<typename T, typename U>

//void assign(ublas::matrix<T>& m, std::size_t r, std::size_t c, U const& data)
//{
//	m.resize(std::max(m.size1(), r + 1), std::max(m.size2(), c +1));//this could be done in a more optimal way
//m(r, c) = data;
//}



/*Calculates U and V factors applying the ACA algorithm.*/
int ACA_new(int m, int n, int* group_o, int* group_t, charge **filaments, double **U, double **V, double tol)
{
	int ooo = 0;

	double lengthdef = sqrt(pow(filaments[0]->fil->length, 2) + pow(filaments[0]->fil->width, 2) + pow(filaments[0]->fil->height, 2));

	int k = 0;
	int JJ = 0;
	int II = 0;
	int *I1 = (int *)calloc(m, sizeof(int));
	int *J1 = (int *)calloc(n, sizeof(int));

	I1[II] = 1;
	double Z1 = 0;
	double E1 = 1;
	double tolerance = tol;
	/*
	if (tol == 0)
	{
	tolerance = 0.00001;
	}
	else if (tol == 1)
	{
	tolerance = 0.0001;
	}
	else if (tol == 2)
	{
	tolerance = 0.001;
	}
	else if (tol == 3)
	{
	tolerance = 0.01;
	}
	else
	{
	tolerance = 0.000001;
	}
	*/

	double *summ2 = (double*)malloc(n * sizeof(double));
	double *summm = (double*)malloc(m * sizeof(double));
	double *error_col = (double*)malloc(m * sizeof(double));
	double *error_row = (double*)malloc(n * sizeof(double));
	if (true)
	{
		while (E1 > Z1*tolerance  && k < n  && k < m)
		{

			lengthdef = filaments[II]->fil->length + filaments[II]->fil->width + filaments[II]->fil->height;

			if (k == 0)
			{
				for (int i = 0; i < n; i++)
				{
					error_row[i] = induct(filaments[group_o[II] - 1]->fil, filaments[group_t[i] - 1]->fil, 0);
				}
				for (int i = 0; i < n; i++)
				{
					if (abs(error_row[i]) > abs(error_row[JJ]))
						JJ = i;
				}
				J1[JJ] = 1;
				for (int i = 0; i < n; i++)
				{
					V[0][i] = error_row[i] / error_row[JJ];
				}
				for (int i = 0; i < m; i++)
				{
					error_col[i] = induct(filaments[group_o[i] - 1]->fil, filaments[group_t[JJ] - 1]->fil, 0);
				}
				for (int i = 0; i < m; i++)
				{
					U[0][i] = error_col[i];
				}
				Z1 = sqrt(pow(frobeniusNormV3(U[k], k + 1, m), 2) *pow(frobeniusNormV3(V[k], k + 1, n), 2));//norm_frobenius(prod(f1, f2));//sqrt(pow(norm_frobenius(U), 2) + pow(norm_frobenius(V), 2)); This needs a seprate function 
				E1 = frobeniusNormV3(U[k], k + 1, m) * frobeniusNormV3(V[k], k + 1, n);
			}
			else
			{
				for (int i = 0; i < m; i++)
				{
					summm[i] = 0;
				}
				for (int i = 0; i < n; i++)
				{
					summ2[i] = 0;
				}
				II = 0;
				while (I1[II] == 1)
				{
					++II;
				}
				for (int i = II; i < m; i++)
				{
					if (abs(error_col[i]) > abs(error_col[II]) && I1[i] == 0)
						II = i;
				}
				I1[II] = 1;


				for (int j = 0; j < k; j++)
				{
					for (int i = 0; i < n; i++)
					{
						summ2[i] += U[j][II] * V[j][i];
					}
				}
				V[k] = (double *)malloc(n * sizeof(double));
				for (int i = 0; i < n; i++)
				{	
					error_row[i] = induct(filaments[group_o[II] - 1]->fil, filaments[group_t[i] - 1]->fil, summ2[i]);
				}
				JJ = 0;
				while (J1[JJ] == 1)
				{
					++JJ;
				}
				for (int i = JJ; i < n; i++)
				{
					if (abs(error_row[i]) > abs(error_row[JJ]) && J1[i] == 0)
						JJ = i;
				}
				J1[JJ] = 1;
				for (int i = 0; i < n; i++)
				{
					if (error_row[JJ] == 0)
						return k;
					else
						V[k][i] = error_row[i] / error_row[JJ];
				}
				for (int i = 0; i < k; i++)
				{
					for (int j = 0; j < m; j++)
					{
						summm[j] += V[i][JJ] * U[i][j];
					}
				}
				U[k] = (double *)malloc(m * sizeof(double));
				for (int i = 0; i < m; i++)
				{
					error_col[i] = induct(filaments[group_o[i] - 1]->fil, filaments[group_t[JJ] - 1]->fil, summm[i]);
				}
				for (int i = 0; i < m; i++)
				{
					U[k][i] = error_col[i];
				}
				double please = multiplyy(k, m, n, U, V);
				Z1 = sqrt(pow(Z1, 2) + please + pow(frobeniusNormV3(U[k], k + 1, m), 2) *pow(frobeniusNormV3(V[k], k + 1, n), 2));
				E1 = frobeniusNormV3(U[k], k + 1, m)*frobeniusNormV3(V[k], k + 1, n);
			} 
			++k;
		}
	}
	return k;
}


